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1.  Summary 


This  report  summarizes  efforts  to  enhance  our  knowledge  of  swirl  injector  dynamics  and  behavior  of 
energetic  additives  that  may  influence  engine  performance  and  stability.  The  swirl  injector  work  was  carried 
over  from  a  prior  contract  and  concluded  about  one  year  ago.  Several  publications  have  stemmed  from  these 
studies  and  we  discovered  fundamental  resonance  conditions  for  these  injectors  that  provides  a  simple  and 
practical  mechanism  to  compute  frequencies  where  substantial  injector/chamber  coupling  may  be  present. 
These  results  were  confirmed  in  experimental  cold  flow  studies  using  a  unique  pulsator  design  developed  in 
our  group. 

The  energetic  particle  studies  have  necessitated  the  development  of  an  entirely  new  model,  based  on  the 
PBE  (Population  Balance  Equation).  The  PBE  has  been  integrated  into  the  GEMS  code  developed  in  Prof. 
Merkle’s  research  group  using  a  fast  Eulerian  method  for  advancing  particle  trajectories  in  time.  The 
DQMOM  (Direct  Quadrature  Methods  of  Moments)  [4]  is  used  for  representation  of  the  particle  size 
evolution  given  an  initial  distribution  of  sizes  in  the  propellants.  Along  with  PBE,  the  harsh  condition  in 
combustion  chamber  and  nozzle  leads  us  to  develop  models  for  laminar  and  turbulent  collision(or 
coalescence)  and  breakup.  The  most  models  which  can  be  found  in  chemistry  and  chemical  engineering 
papers  are  limited  only  by  the  turbulent  viscous  effect  in  a  low  turbulence  case.  Unlike  these  investigations, 
the  drops  in  a  combustion  chamber  are  exposed  to  a  highly  turbulent  flow  and  consequently  inertial  effects 
of  the  drops  induced  by  larger  density  of  particles  than  the  surrounding  gas  are  important.  So,  we  divide  the 
collision  and  breakup  mechanisms  into  four  regimes  and  each  regimes  are  modeled:  laminar  hydrodynamic 
collision/breakup,  laminar  aerodynamic  collision/breakup,  turbulent  hydrodynamic  collision/breakup,  and 
turbulent  aerodynamic  collision/breakup.  Here,  the  term,  hydrodynamic,  means  the  shearing  motion  of 
surrounding  fluid  is  the  main  source  for  collision  and  breakup  and  the  term,  aerodynamic,  means  the 
velocity  difference  between  the  particle  and  the  surrounding  fluid  is  the  main  source  for  collision  and 
breakup. 
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Technology  Transfer 

The  research  group  at  Purdue  is  supporting  a  variety  of  developments  throughout  the  industry.  Under 
NASA  sponsorship,  models  that  were  initially  created  in  the  AFOSR  program  are  being  used  to  assess  the 
forced  response  of  plain  orifice  “pressure  atomizers”  under  a  wide  range  of  conditions.  The  models  are 
being  incorporated  into  the  industry-standard  Rocket  Combustor  Interactive  Design  (ROCCID)  code  that  is 
used  by  NASA  MSFC,  U.S.  Air  Force,  and  numerous  propulsion  contractors.  The  models  will  substantially 
improve  the  basic  treatment  of  these  atomizers  and  the  Purdue  team  is  working  closely  with  Sierra 
Engineering  on  implementation  of  the  new  models.  Results  from  current  AFOSR-sponsored  efforts  in  the 
dynamics  of  swirl  injectors  has  also  been  transmitted  to  NASA  officials  as  well  as  prior  simulations  of  shear 
coaxial  injectors  that  are  of  great  interest  for  new  Crew  Exploration  Vehicle  propulsion.  Our  team  works 
closely  with  small  companies  including  Sierra  Engineering  and  INSpace  LLC  to  provide  recommendations 
on  injector  designs.  We  have  also  provided  inputs  on  gas/gas  injectors  for  potential  application  to  lunar 
transfer  vehicles  under  sponsorship  from  entities  affiliated  with  Kistler  Aerospace.  Currently,  we  are 
working  on  nonlinear  dynamics  of  swirl  injectors  and  hope  to  be  able  to  create  a  submodel  for  ROCCID  that 
would  incorporate  these  results  and  permit  the  code  to  assess  a  whole  new  class  of  injectors.  A  comparable 
submodel  for  shear  and  swirl  coaxial  injectors  is  also  under  development  under  NASA  sponsorship, 
although  this  is  a  rather  low-level  effort  at  present. 


2.  Research  Objectives 


The  understanding  of  the  complex  combustion  phenomena  present  in  liquid  rocket  engines  begins  with 
the  fundamental  process  of  fuel  and  oxidizer  jet  atomization.  The  objective  of  this  research  has  been  to 
develop  a  model  to  account  for  agglomeration  and  breakup  of  oxides  of  energetic  particles  used  in  solid  and 
liquid  propellants.  The  focus  of  the  work  is  to  understand  the  role  these  particles/droplets  play  in 
combustion  stability  and  overall  engine  performance. 

At  present  there  is  substantial  interest  in  studying  the  effect  of  nanoenergetic  ingredients  that  have 
recently  become  available  as  a  result  of  manufacturing  process  advances.  The  nanomaterials  provide  for 
dramatic  reductions  in  the  time  required  for  combustion  events,  but  pose  challenges  in  terms  of 
incorporation  in  the  propellant  and  in  terms  of  the  tremendously  dense  cloud  of  particles  formed  as  a  result 
of  the  combustion  process.  In  contrast  to  micron-sized  particles  used  in  prior  research,  agglomeration  events 
become  much  more  prevalent  when  nanomaterials  are  employed.  Breakup  of  agglomerates  and  potential 
wall  impingement  in  the  throat/nozzle  region  are  also  important  processes  that  effect  overall  motor 
performance  via  the  two-phase  flow  loss  and  nozzle  erosion.  For  this  reason,  a  computational  tool  has  been 
developed  to  assess  agglomeration  and  breakup  to  assess  particle  size  development  for  arbitrary  combustor 
type  (solid  motor,  liquid  or  hybrid  rocket).  The  objectives  of  the  study  is  to  provide  a  theoretical  basis  for 
performance  changes  attributed  to  use  of  the  nanoenergetic  materials  and  to  assess  design  parameters  that 
influence  such  performance. 


3.  Status  of  Research  -  Development  of  Energetic  Particle 
Agglomeration  and  Breakup  Methodology 


I.  Introduction 


Metal  additives  are  used  to  enhance  the  energy  of  solid  propellants  and  are  also  given  consideration  for  loading  in 
liquid  slurries  for  the  same  reason.  Historically,  micron-sized  particles  have  been  used  for  this  application,  but  the 
recent  large-scale  manufacture  of  nano-sized  particles  by  numerous  vendors  changes  the  overall  number  of  particles 
and  their  average  spacing  dramatically.  For  this  reason,  collision/agglomeration  processes  that  may  have  been 
neglected  in  the  past  could  be  of  significant  importance  in  these  flows.  In  case  of  aluminum  loaded  solid  propellants, 
Gany  et  al.  [1]  have  experimentally  investigate  the  A1/A1203  agglomerates  forming  on  the  propellant  surface  and 
observed  a  mean  of  about  250  pm  when  a  propellant  contains  6  pm  aluminum  particles.  More  recently,  Najjar  et  al.  [2] 
have  referred  to  Sabnis  et  al.  [3]  and  indicated  that  the  typical  values  of  the  drop  size  distribution  entering  the 
combustion  chamber  are  a  mean  of  150  pm  for  larger  Aluminum  particle  and  1.5  pm  for  smaller  Aluminum  Oxide 
particle,  which  is  bimodal.  Nano-aluminum  loaded  solid  propellants  exhibit  significantly  different  agglomeration  near 
the  surface  and  thus  the  size  of  agglomerates  leaving  the  burning  surface  can  be  significantly  smaller.  Much  less  is 
known  about  the  potential  for  agglomeration  of  particle  loaded  liquid  slurries/gels,  but  regions  of  high  shear  that  are 
present  due  to  mixing  processes  could  presumable  provide  substantial  opportunities  for  agglomeration  to  occur.  The 
idealized  concept  of  particle  size  variation  in  rocket  chamber  is  illustrated  in  Fig.  1.1. 


Fig.  1.1  Illustration  of  the  simple  concept  of  droplet  size  variation  in  rocket  chamber 


A  simple  analysis  of  the  particle  to  particle  distance  in  rocket  chamber  can  show  the  substantial  opportunities  for 
agglomeration  due  to  collision  of  particles.  As  Najjar  et  al.  [2]  have  indicated  that  the  burning  of  20%  aluminum 
loaded  propellant  of  solid  rocket  booster  results  in  approximately  1015  droplets  of  a  mean  diameter  of  100  pm  in  a  core 
volume  of  63  m3.  Following  Friedlander  [4],  the  average  center-to-center  distance  between  two  adjacent  particles 
distributed  randomly  is  given  by  0.55396A/oo"13  where  TV is  a  number  density  of  particles.  Given  the  number  density 
using  the  data  of  Najjar  et  al.  (Aoo=number  of  droplets/contained  volume),  the  average  distance  between  two  adjacent 
particles  is  approximately  22  pm.  Considering  the  droplets  of  a  mean  of  100  pm,  and  that  the  overall  distance  traveled 
is  of  the  order  of  10’s  of  meters,  it  is  inevitable  that  collisions  will  occur. 

The  drops  entering  the  chamber  have  a  substantial  effect  on  the  rocket  motor  efficiency.  The  increase  of  the 
specific  impulse  and  damping  of  chamber  instability  may  be  desirable  effects.  However,  slag  accumulation,  nozzle 
erosion,  and  significant  exhaust  signature  are  disadvantages  of  aluminum  loaded  propellant.  The  particle  phase 
characteristics,  especially  the  number  density  (or  mass  concentration)  and  the  drop  size  may  be  thought  as  governing 
parameters  in  assessing  these  effects.  As  particles  exit  the  nozzle  with  velocities  less  than  the  gas  depending  on  their 
size  and  drag  characteristics  a  two-phase  flow  loss  is  always  present  in  gas/particle  nozzle  flows.  Therefore,  the 
prediction  of  particle  phase  characteristics  is  of  high  importance  in  quantifying  two-phase  flow  losses  and  ascertaining 
performance  advantages. 

The  past  and  present  studies  of  two  phase  flow  inside  the  rocket  chamber  have  focused  on  the  effects  of  the 
droplet  on  the  gas  flow  by  two-way  coupling  [2,  5,  6]  and  the  effects  of  the  gas  flow  on  the  particle  phase  by  one-way 
coupling  [6].  However,  none  of  these  studies  have  been  focusing  on  the  effects  of  the  collision  and  breakup  of  the 
droplets  and  consequent  drop  size  change.  Although  Najjar  et  al.  [2]  have  included  the  collision  effects  in  assessing  the 
drop  mass  change,  the  collision  efficiency  in  their  model  is  simply  set  as  a  constant,  0.25. 


The  flow  in  a  large  rocket  chamber  can  experience  highly  shearing  motion  due  to  its  mean  value  change  and 
highly  turbulent  motion  at  the  same  time.  The  high  Reynolds  number  and  the  complex  geometry  of  solid  rocket 
chamber  leads  to  the  locally  complex  flow  motion  and  two  adjacent  particles  can  be  easily  intercepted  by  the  turbulent 
motion  of  flow.  In  addition,  the  highly  shearing  motion  of  mean  flow  near  boundary  layer  can  result  in 
collision/breakup.  Thereby,  stochastic  collision  and  breakup  events  can  be  one  of  the  governing  mechanism  of  the 
particle  to  particle  interaction  in  a  rocket  chamber  and  collision  and  breakup  due  to  mean  flow  motion  can  be  another 
governing  mechanism. 

The  coalescence  and  breakup  process  of  drops  [7,  8]  and  bubbles  [9]  have  been  investigated  in  the  chemistry  and 
chemical  engineering  communities.  The  modeling  of  the  coalescence  and  breakup  processes  in  an  agitated  vessel  have 
been  important  topic  in  chemistry  to  assess  the  mixing  effects.  Their  interests  are  usually  limited  only  by  the  turbulent 
viscous  effect  in  a  low  turbulence  case.  Unlike  the  two  immisicible  liquids  in  an  agitated  vessel,  the  drops  in  a 
combustion  chamber  are  exposed  to  a  highly  turbulent  flow  and  consequently  inertial  effects  of  the  drops  induced  by 
larger  density  of  particles  than  the  surrounding  gas  are  important.  These  factors  lead  to  difficulties  in  using 
coalescence  and  breakup  models  developed  in  chemistry  but  these  models  can  be  a  good  starting  point  in  current 
modeling. 

Thus,  the  stochastic  collision  and  breakup  are  addressed  here  and  the  collision/breakup  in  laminar  flow  and 
combination  of  the  mean  flow  effects  and  turbulent  flow  effects  are  discussed  too.  We  divide  the  collision  and  breakup 
mechanisms  into  four  regimes  and  each  regimes  are  modeled:  laminar  hydrodynamic  collision/breakup,  laminar 
aerodynamic  collision/breakup,  turbulent  hydrodynamic  collision/breakup,  and  turbulent  aerodynamic 
collision/breakup.  Here,  the  term,  hydrodynamic,  means  the  shearing  motion  of  surrounding  fluid  is  the  main  source 
for  collision  and  breakup  and  the  term,  aerodynamic,  means  the  velocity  difference  between  the  particle  and  the 
surrounding  fluid  is  the  main  source  for  collision  and  breakup.  More  details  on  each  term  are  given  in  Chapter  2. 

Along  with  the  collision/breakup  models,  to  assess  the  particle  phase  velocity  field  while  holding  the  reasonable 
computational  efficiency,  an  Eulerian-Fast  (or  Equilibrium)  Eulerian  two-phase  methodology  is  chosen  and  the  direct 
quadrature  method  moment  (DQMOM)  approximation  is  applied  to  the  population  balance  equation  (PBE)  is  used  to 
model  the  coalescence  and  breakup.  The  details  of  methodology  are  provided  in  Chapter  2. 

The  objective  of  the  current  study  is  to  develop  models  for  the  collision  and  breakup  processes  applicable  to  a 
simulation  of  the  two  phase  flow  in  a  rocket  chamber  and  carry  a  test  simulation  in  a  typical  rocket  chamber  and 
attached  converging-diverging  nozzle.  For  this  purpose,  Computations  were  performed  on  a  typical  converging- 
diverging  nozzle  attached  to  a  rocket  motor.  The  MMD  (Mass  Mean  Diameter)  was  predicted  according  to  different 
droplet  characteristics  and  pressure  at  nozzle  inlet  and  the  scales  of  nozzle.  To  validate  the  models,  the  predicted 
results  are  compared  to  Hermsen  [10] ’s  empirical  correlation  which  predicts  the  particle  size  at  the  exit  plane  of  SRM 
nozzle.  The  results  are  reasonably  agreed  with  the  empirical  correlation.  However,  the  simulation  is  very  sensitive 
with  the  initial  droplet  condition  (i.e.  mean  diameter  and  standard  deviation  of  number  distribution),  therefore,  the 
initial  conditions  of  droplets  should  be  chosen  very  carefully. 


II.  Physical  Modeling 


2.1  Flow  field  description  -  Navier-Stokes  equation 


The  2-D  unsteady  Navier-Stokes  equations  for  the  Newtonian  viscous  carrier  fluid  are  applicable  under  the 
continuum  condition.  The  flow  field  is  described  by  mass,  momentum  and  energy  conservation  laws  complemented  by 
an  appropriate  equation  of  state  and  additional  constitutive  relations.  Two  turbulence  equation  from  the  k-co  model  of 
Wilcox  [11]  are  added  to  the  conservation  form  of  the  Navier-Stokes  equations  without  any  body  forces  and  source 
terms  induce  by  the  particle  phase: 
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where  the  vectors,  Q  ,  E  ,  V  ,  and  H  are  given  by 
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The  Xi  and  ut  represent  Cartesian  coordinates  and  velocity  components,  P  and  p  represent  the  pressure  and  density,  h° 
is  the  stagnation  enthalpy,  and  K  is  the  thermal  conductivity.  The  A,  A  ,  B,  and  5*  are  closure  constants  for  Wilcox 
turbulence  model,  Ty  is  the  Reynolds  stress  tensor,  and  v ^  is  the  turbulence  eddy  viscosity.  These  terms  are  given  as 
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The  gas  phase  governing  equations  is  solved  under  the  platform  of  Generalized  Equation  and  Mesh  Solver  (GEMS) 
code  [12],  which  uses  contemporary  numerical  methods  to  solve  coupled  systems  of  partial  differential  equations. 


2.2  Particle  phase  modeling  -  Population  Balance  Equation  (PBE) 


—  +  V-(upn)  =  V-(DsVn)  +  S 


(2.4) 


2.2.1  QMOM  (Quadrature  Methods  of  Moments) 

The  advection-diffiision  equation  for  the  number  density  field  is  given  by 

dn 
dt 

where  n  is  is  the  particle  number  density,  Ds  is  the  diffusion  coefficient,  and  S  is  the  source  term  corresponding  to 
coagulation  and  breakup.  In  a  high  Reynolds  number  or  shearing  flow,  the  diffusion  term  can  be  ignored  and  the 
advection-diffusion  equation  becomes  a  form  similar  to  Smoulchowski’s  equation  [13]  which  is  usually  referred  as  the 
population  balance  equation. 

Using  a  one-way  coupling  approach,  no  mass,  momentum,  and  energy  interchange  is  considered.  The  particle 
phase  is  also  assumed  to  be  in  thermally  equilibrium  state.  The  equation  constructing  agglomeration/breakage  models 
for  the  dispersed  phase  is  the  population  balance  equation  for  the  particle  number  density  which  is  as  follows  [14]: 
dn(v,t) 
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This  equation  assumes  the  same  velocity  of  particle  phase  as  the  surrounding  fluid.  This  equation  expresses  the  fact 
that  there  is  break-up  and  coalescence  of  droplets  in  the  flow  in  the  absence  of  interactions  with  walls.  The  term 
n(y,t)  is  the  number  density  function  of  the  particle  volume  v  and  u  is  the  velocity  of  the  carrier  fluid  velocity  due  to 
Stokesian  particle  assumption.  Here,  a,  is  the  collision  efficiency  between  particles  with  volume  v  and  v  ,  a2  is  the 

coalescence  efficiency,  (3  is  the  volume  based  collsion  kernel  that  describes  the  frequency  that  particles  of  volume  v 

* 

and  v  collide,  a  is  the  fragment  distribution  function,  and  b  is  the  volume -based  breakage  kernel  that  is  the 
frequency  of  breakage  of  a  particle  of  volume  v  [14].  The  first  term  on  the  right-hand  side  represents  the  formation  of 
volume  v  by  collision  and  the  second  term  represents  the  loss  of  the  volume  v  by  collision.  The  third  term  represents 
the  formation  of  volume  v  by  break-up  and  the  last  term  represents  the  loss  of  volume  v  by  break-up. 

Solving  this  equation  directly  will  require  large  computational  power  due  to  the  presence  of  a  large  number  of 
classes  of  particles.  In  addition,  the  source  terms  in  the  equation  represent  that  the  equations  for  each  phase  are  highly 
coupled  by  each  other.  Therefore,  the  simplification  of  the  governing  equations  is  highly  required.  This  can  be 
achieved  by  QMOM  (Quadrature  Method  of  Moments)  developed  by  Mcgraw  [15]  which  is  a  powerful  technique  to 
determine  the  evolution  of  the  lower-order  moments  of  the  distribution  by  a  quadrature -based  approximation.  Wang  et 
al.  [14]  have  successfully  applied  this  approach  in  Taylor  Coutte  flow,  and  Marchisio  et  al.  [16]  have  showed  that  this 
approach  leads  to  very  small  error  comparing  to  discretized  population  balance  equation  (DPB).  Wang  et  al.  [14]’s 
length  based  QMOM  approximation  process  of  PBE  is  summarized  here.  The  QMOM  starts  from  defining  the 
moments  and  taking  quadrature  approximation  as  follows: 
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where  l3  —  v,  z  —  3k 


The  term  N  is  the  order  of  the  quadrature  formula  and  v  is  the  particle  volume.  Accordingly,  m0  is  the  total  particle 
number  density  and  M\  is  the  total  particle  volume  concentration  (same  as  particle  volume  fraction).  Applying  the 
length-based  definition  of  moments  to  the  transport  equation  of  the  particle  density  gives  (superscript  '  is  omitted 
here): 
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The  collision  efficiency  is  a*.  =  cq(/.,/y.) ,  the  coalescence  efficiency  is  a*-  =  a2 (f , T ) ,  the  collision  frequency  of 
drops  of  length  lt  and  l,  is  (3y  =  /?(/.,  /  .) ,  the  breakage  frequency  of  drops  of  length  //  is  bi  =  b(l) ,  and  the  daughter 
drop  probability  density  function  for  the  binary  fragmentation  is  given  by  [14] 

a  =2(3-z)/3/z 


(2.8) 


The  weights,  wt  and  abscissas,  /•  are  found  via  using  of  the  product-difference  (PD)  algorithm. 

The  product  difference  (PD)  algorithm,  which  is  used  to  find  weights  ( w. )  and  abscissas  ( /• )  from  the  moments 

( mz )  while  solving  PBE,  is  given  by  Mcgraw  [15]  and  Wang  et  al.  [14]  and  it  is  summarized  here. 

The  first  step  is  to  obtain  a  matrix  P  as  follows: 
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From  the  obtained  matrix  P,  the  coefficients  (a, )  are  given  as  follows: 

=0 

=  Px,M  1  (PuPu-i )  for  i=2> . ’ 27V 

Then,  a  symmetric  tridiagonal  matrix  S  is  obtained  with  the  following  diagonal  ( sd  i )  and  co-diagonal  ( scd  i ) 
components: 


Sd,i  ~  a2i  + 


for  /=  1, . ,N 

for  i=l, . N- 1 


(2-11) 


After  the  symmetric  tridiagonal  matrix  is  obtained,  the  weights  and  abscissas  are  obtained  by  finding  its  eigenvalues 
and  eigenvectors.  The  eigenvalues  of  the  matrix  S  are  the  abscissas  and  the  weights  are  given  by 
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where  viX  is  the  first  component  of  eigenvector  za  .  The  eigenvalues  and  eigenvectors  are  found  by  QL  algorithm  [17]. 


2.2.2  DQMOM  and  a  fast  Eulerian  approach 

The  number  of  drops  in  a  chamber  is  typically  very  dense  as  explained  above  and  then  the  particle  phase  is  treated 
as  continua.  Therefore,  the  particle  phase  can  be  described  via  an  Eulerian  approach.  For  more  simplicity  and 
numerical  efficiency,  a  fast  (or  equilibrium)  Eulerian  approach  [2,  18]  is  used,  such  that  mass  and  momentum 
conservation  are  automatically  satisfied.  In  this  approach,  the  particle  phase  velocity  is  handled  as  a  field  variable 
which  is  given  by 

UP=U-T—  (2.13) 

Dt 


where  U p  and  U  are  the  particle  phase  and  gaseous  phase  velocity  vectors,  respectively.  The  term  t  is  the  relaxation 
time  of  the  particle  and  D I  Dt  is  the  material  derivative  in  the  Eulerian  view. 

Because  the  particles  have  larger  density  than  the  gaseous  phase,  a  distribution  over  particle  velocities  is  needed  to 
be  considered.  A  multivariate  number  distribution  function  n  depends  on  /,  Uh  xh  t  which  can  be  denoted  as  n(l,  Uh  xh 
t).  In  this  case,  the  transport  equation  is  given  by 
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which  was  proposed  by  Williams  [19]  for  the  LHS. 

To  reduce  the  number  of  variables,  the  averaged  number  distribution  function  and  the  averaged  drop  phase 
velocities  can  be  given  as  follows: 
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Integrating  the  equation  (2.10)  over  the  velocity  assuming  a  Dirac  delta  function  of  velocity  distribution  yields  the 
following  PBE: 


n(l,xi,t)  +  —(Upin(l,xi,t))  =  S(l,xi,t) 


dt 


(2.16) 


which  is  same  as  the  equation  (2.4).  Following  Marchisio  and  Fox  [20],  the  particle  size  distribution  function  can  be 
treated  as  a  sum  of  Dirac  delta  functions: 
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Substituting  the  equation  (2.13)  into  the  equation  (2.12)  and  integrating  the  equation  (2.12)  gives  the  following  the 
PBE  approximated  by  DQMOM: 
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where  ^  is  the  weighted  abscissas  defined  by  £  =  w  l  .  The  source  terms  can  be  obtained  by  solving  the  following 


equation. 
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where  Sm,k  is  the  source  term  obtained  in  QMOM  case  and  it  is  given  in  the  RHS  of  the  equation.  The  source  terms, 
S  and  S.  t/ ,  can  be  obtained  from  the  linear  system  Ax=b  which  each  matrix  is  defined  as  follows: 
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2.3  Stochastic  coalescence/breakup  modeling 


2.3.1  Collision  frequency  kernel 

2. 3. 1.1  Spherical  formulation  of  collision  kernel 

Saffman  and  Turner  [21]  have  studied  the  collision  frequency  kernel  and  presented  two  formulations  of  the 
collision  kernel  which  are  spherical  formulation  and  cylindrical  formulation.  Comparing  these  two  formulations,  the 
collision  kernel  in  cylindrical  formulation  is  the  cylindrical  volume  passing  through  the  effective  collision  circle  per 
unit  time  (in  other  words,  the  fluid  volume  flux  across  the  effective  collision  area)  and  the  collision  kernel  in  spherical 
formulation  is  the  volume  of  fluid  across  the  collision  sphere  surface  (volume  flux  across  the  collision  sphere  surface). 
The  cylindrical  formulation  is  possible  in  the  special  case  of  uniform  shear  flow  which  is  same  case  with 
Smoulchowski  [13]  and  the  more  general  way  will  be  the  spherical  formulation  because  the  relative  velocity  between 
particles  depends  on  the  orientation  of  the  collision  radius  Rc  as  it  is  described  in  Wang  et  al.  [22]. 


y 


Fig. 2.1  Schematic  of  collision  of  two  droplets  of  radii  r7  and  r2  in 
Spherical  formulation;  the  collision  radius  Rc  is  the  sum  of 
radii  r7  and  r2,  the  relative  motion  follows  the  streamlines 

The  Saffman  and  Turner  [21] ’s  spherical  formulation  is  described  in  Figure  2.1.  Considering  two  particles  of  radii 
rx  and  r2,  the  moving  particle  is  the  particle  of  radius  r2  supposing  the  particle  of  radius  rx  as  a  fixed  central  particle. 
Assuming  there  is  no  distortion  of  flow  field  due  to  the  existence  of  the  particle,  the  particle  r2  are  moving  along  the 
streamlines.  Defining  the  collision  sphere  as  a  sphere  of  radius  Rc=rx+r2  centered  on  the  fixed  central  particle,  the 
collision  frequency  of  the  fixed  central  particle  is  the  flux  of  the  fluid  having  the  velocity  which  is  same  as  the  relative 
velocity  between  two  particles,  multiplied  by  the  number  density  of  the  moving  particles.  This  flux  should  be  induced 
by  the  relative  velocity  which  is  inwardly  normal  to  the  collision  sphere  because  this  directional  component  of 
velocity  is  only  the  component  causing  the  collision.  Denoting  the  unit  vector  outwardly  normal  to  the  collision  sphere 

(radial  direction  of  the  collision  sphere)  as  nr  and  the  relative  velocity  inwardly  normal  to  the  collision  sphere  as  W  , 
the  flux  Ji  across  the  collision  sphere  is  given  by 

J,=-Jw-nrdA  (2.21) 

The  negative  sign  is  given  because  the  dot  product  between  this  velocity  vector  and  outward  normal  vector  is  negative. 
Supposing  that  the  particles  are  distributed  in  the  flow,  the  collision  frequency  Nc  which  is  the  total  number  of 
collision  between  particles  of  the  number  densities  nx  and  n2  in  unit  volume  and  unit  time  is  given  by 

Nc=-nxn2jw-nrdA  (2.22) 

Thus,  the  collision  frequency  function  (or  collision  kernel)  Pi  for  the  laminar  flow  is  given  by 

fy  —  J,  ——Jw  nr  dA  (2.23) 

where  dA  is  the  area  element  on  the  surface  of  a  sphere. 

Developing  further  for  the  turbulent  flow,  when  the  particles  are  randomly  distributed  and  their  fluctuating  radial 
velocity  component  is  wr  (sign  of  this  component  are  not  decided  yet,  the  effect  of  its  mean  component  <Wr>  is  not 
considered),  the  mean  flux  J,  towards  the  collision  sphere  is  given  by 


(2.24) 
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Assuming  the  outward  and  inward  fluxes  across  the  collision  sphere  are  equal,  this  assumption  can  be  expressed  by 

—  j  wr  dA—  J  wr  dA  (2.25) 


wr  <0  wr  >0 

Thus,  the  mean  flux  J,  towards  the  collision  sphere  is  given  by 
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For  isotropic  turbulence,  the  collision  kernel  /3t  is  given  by 

A  =2^c(kl) 

which  is  the  same  form  as  Saffman  and  Turner  [21]’s. 


(2.26) 


(2.27) 


2. 3. 1.2  Collision  frequency  kernel  -  Turbulent  flow 

The  present  study  is  focused  on  the  collision  (or  coalescence)  and  breakup  of  two  unequal  spherical  drops. 
Concerning  the  hydrodynamics  between  two  drops,  the  assumption  of  two  equal  drops  will  make  it  easy  to  analyze 
flow  motion.  The  collision  of  equally  assumed  two  drops  is  described  in  Figure  2.2.  Following  Chesters  [23],  two 
unequal  spherical  drops  can  be  characterized  by  two  equal  drops  of  equivalent  radius,  Req,  which  is  given  by 

2  RR 
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where  Rt  and  Rj  are  radii  of  two  unequal  drops.  When  a  drop  of  Req  is  smaller  than  Kolmogorov’s  length  scale,  the 
drop  is  considered  as  it  is  in  the  viscous  subrange  of  turbulence  and  a  drop  larger  than  Kolmogorov’s  length  scale  is 
considered  as  it  is  in  the  inertial  subrange.  The  Kolmogrov’s  length  scale  rj  is  given  by 


(2.29) 


where  v  is  the  kinematic  viscosity  of  the  fluid. 


Fig. 2. 2  Schematic  of  droplet  collision;  left  is  for  two  unequal  droplets  collision, 
right  is  for  two  equally  assumed  droplets  collision  (Ri=R^=Req) 


2. 3. 1.3  Relative  velocity  between  two  drops 

Before  discussing  the  collision  frequency  function,  it  is  convenient  to  consider  the  relative  velocity  between  two 
colliding  drops.  In  analogy  with  Williams  and  Crane  [24],  the  relative  velocity  between  two  particles  can  be  thought 
as  it  is  induced  by  two  major  effects:  the  effect  of  velocity  gradient  of  the  carrier  fluid  between  two  particles  and  the 
effect  of  different  inertial  response  of  particles  of  different  radius  to  the  movement  of  carrier  fluid.  According  to  these 
considerations,  the  relative  velocity  between  two  particles  can  be  constituted  by  the  effects  of  velocity  gradient  and 

different  inertial  response.  It  is  supposed  that  two  droplets  within  the  fluid  have  velocities  UP, i  and  U p,i  before  they 

collide.  The  carrier  fluid  surrounding  these  drops  have  velocities  U i  and  U 2 .  When  the  slip  velocities  between  the 

particle  and  the  carrier  fluid  are  denoted  by  Qx  —  U p,  1  —  U 1  and  Q2  =  U p, 2  — Ui  ,  the  relative  velocity  vector 

W  =  Up, 2  —  U p,\  can  be  expressed: 


W  =  Wi+Ws 


(2.30) 


=  +  U2-U1 

The  first  term  on  RHS  Wi  represents  the  amount  of  velocity  difference  induced  by  inertial  effects  of  large  density 

particles  and  the  second  term  Ws  represents  the  velocity  difference  induced  by  the  velocity  gradient  (or  strain  rate)  of 
the  carrier  fluid.  The  modeling  of  each  term  in  the  laminar  flow  will  be  accomplished  in  the  future. 

Considering  the  collision  of  two  particles  of  radius  rt  and  rp  the  responsible  component  of  relative  velocity  to  the 
collision  is  only  the  component  in  the  direction  of  the  centerline  which  connects  the  center  of  two  particles.  Figure  2.3 
illustrates  which  components  of  the  relative  velocity  are  related  to  the  collision.  Supposing  the  fixed  central  particle  of 

radius  r7,  the  moving  particle  rt  has  velocity  Wi  and  Ws  of  which  Wi,r  and  W s,r  makes  two  particles  to  approach 

each  other  whereas  Wi,t  and  Ws,t  cause  movement  away  from  the  fixed  central  particle.  Therefore,  only  the  velocity 
components  that  induce  approaching  motion  must  be  considered  and  the  velocity  components  causing  movement 
away  must  be  neglected  from  the  consideration. 


Fig.2.3  Schematic  of  components  of  relative  velocity  between 
two  particles  of  radii  rt  and  r7 


Considering  the  turbulent  flow,  the  mean  square  of  the  relative  velocity  is  given  by  [21] 

iw-w)  =  (wi  -Wi\  +  l Ws 


(2.31) 

It  is  assumed  that  ws  is  statistically  independent  of  the  slip  velocities  qx  and  q2  .  Using  the  notation 
ws  —  ws  x ,  Ws  y ,  ws  z  9  the  second  term  on  the  RHS  is  given  by 

(ws  ■  ws^l  =  ( w2Sx )  +  (w2Sy )  +  ( wlz ) 

Using  the  mean  square  of  the  velocity  gradient  in  viscous  subrange  which  is  given  by 
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For  isotropic  turbulence,  it  has  been  shown  that  [26] 
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(2.32) 

(2.33) 


(2.34) 


(2.35) 


(2.36) 


This  form  is  same  as  the  velocity  gradient  term  in  the  analysis  of  Saffman  and  Turner  [21]  considering  the  inertial 
effect  together  with  the  effect  of  velocity  gradient  term. 

To  evaluate  the  first  inertial  effect  term  (wi  •  wi  \ ,  it  is  assumed  that  the  velocities  of  the  carrier  fluid  near  two 


adjacent  particles  are  the  same.  This  assumption  is  valid  when  the  two  particles  are  smaller  than  the  smallest  eddy,  in 
other  words,  it  is  valid  when  they  are  in  viscous  subrange.  This  assumption  is  equivalent  to  the  assumption  of  Saffman 
and  Turner  [21]  which  is  used  in  evaluating  the  same  term  as  the  first  inertial  effect  term  of  this  paper.  This 

assumption  results  in  the  assumption  that  the  correlation  coefficient  between  qx  and  q2  is  unity  where  qx  and  q2  are 

the  fluctuating  parts  of  the  relative  velocity  between  particle  and  its  surrounding  fluid  of  particle  1  and  2  respectively. 
The  evaluation  of  this  term  will  be  discussed  in  the  future  work.  Using  the  assumption,  the  first  inertial  effect  term  can 
be  given  by 

/wi-Wl\  =  (  Up,2—Up,\  •  Up,2—Up,\  \  (2.37) 


li  p  2  U  p  1  •  U  p  2  u  n  1 


Therefore,  the  i-direction  mean  square  relative  velocity  between  two  particles,  (Wj .  j  can  be  expressed: 

K  >  Ru )  +  «2,  >  -  2  {upXiupXi )  (2.38) 

This  equation  is  the  starting  point  of  Williams  and  Crane  [24]’ s  analysis  of  the  fluctuating  relative  motion  of  two 
particles  induced  by  slip  motion  between  the  particle  and  fluid.  Saffman  and  Turner  [21]  also  have  derived  the  mean 
square  of  the  relative  velocity  (w-w)  considering  the  effects  of  velocity  gradient  and  the  inertial  effect  in  viscous 

subrange.  However,  in  fundamental  assumption  of  their  approach,  the  relaxation  time  is  smaller  than  the  time  scale  of 
smallest  eddy.  So,  the  term  for  the  inertial  effect  cannot  be  used  in  our  case.  Instead,  Williams  and  Crane  [24]  have 
derived  the  term  of  inertial  effect  for  the  small  drops  which  is  not  limited  by  the  small  relaxation  time.  In  Williams 
and  Crane  [24]’ s  analysis,  the  particle  motion  is  described  by  the  simplified  Tchen  [27]’ s  force  balance  equation 
ignoring  the  added  mass,  Basset  history,  and  gravitational  acceleration  terms: 

^  d±+U-L  =  -  (239) 

dt  t  2  pp  +  pg  dt  t  9  ji 

where  ut  and  upJ  are  the  fluctuating  parts  of  the  fluid  and  particle  velocities  in  i-direction,  pg  and  pp  are  the  densities  of 
the  fluid  and  particle,  and  r  is  the  relaxation  time  of  the  particle  of  radius  r.  Stoke ’s  drag  law  is  applied  here 
assuming  the  size  of  the  particle  is  small  enough.  For  pp  »  pg, ,  the  first  term  on  RHS  can  be  neglected.  In  analogy  with 
Levins  and  Glastonbury  [28],  using  a  more  accurate  form  of  the  wave -number  spectrum  which  includes  the  influence 
of  turbulence  within  the  viscous  subrange,  the  mean  square  of  the  particle  velocity  is  obtained  by  Williams  and  Crane 
as  follows: 


(2.39) 
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The  mean  square  fluctuating  velocity  (uj  is  given  in  terms  of  turbulent  kinetic  energy  k  assuming  isotropic 

turbulence.  The  8  is  the  non-dimensional  relaxation  time  of  drop  of  radius  r,  z  is  the  particle  relaxation  time,  and  Lf  is 
the  longitudinal  integral  length  scale  which  is  approximated  by  0.5 L  where  L  is  the  integral  length  scale  of  the  largest 
energy-containing  eddy  approximated  by  kxl2l{B*co)  in  k-co  model.  Finally,  X  is  the  Taylor’s  microscale  length. 

Starting  from  Panchev  [29] ’s  integrated  form  of  Tchen’s  equation  of  motion,  Williams  and  Crane  derived  the 
covariance  (up ,  tup  2 .  ^  using  the  more  accurate  wavenumber  spectrum  for  the  small  particles  which  satisfies  the 
condition  6  «  1  (which  imposes  that  the  particle  is  in  viscous  subrange).  Thus,  using  the  mean  square  particle 
velocity  and  the  covariance  terms,  the  i-direction  mean  square  relative  velocity  between  two  particles  (wr)  for 
viscous  subrange  can  be  expressed: 

For  large  particles  which  satisfy  the  condition  9  »  1  (which  imposes  the  particle  is  in  energy  containing  region),  it 
was  shown  that: 


(2.41) 
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(2.42) 


They  have  successfully  showed  that  the  mean  square  relative  velocity  for  viscous  subrange  goes  to  Saffman  and 
Turner  [21]’s  inertial  term  at  lower-limit  of  small  relaxation  time  and  the  mean  square  relative  velocity  for  energy 
containing  region  goes  to  Abrahamson  [30] ’s  term  at  higher-limit  of  large  particle  size.  They  also  derived  the 
universal  solution  of  the  mean  square  relative  velocity  which  can  be  used  in  inertial  subrange,  it  is  given  by: 

,  I  1  +  0,  +  0. 

(d.+O.f- 46.9.  - ^ - J— 

K>  -«)— _ (2,43) 
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The  universal  solution  of  Williams  and  Crane  approaches  Abrahmson  [30] ’s  mean  square  relative  velocity  for  energy- 
containing  range  at  higher-limit.  It  should  be  noted  here  that  their  universal  solution  does  not  approach  Saffman  and 
Turner’s  inertial  term  as  it  is  explained  in  Kuris  and  Kusters  [31].  However,  the  divergence  for  very  small  particle  is 
appreciable  as  it  is  shown  by  Williams  and  Crane  [24]. 

Finally,  the  mean  square  relative  velocity  induced  by  different  inertial  response,  (wi-wi\9  can  be  calculated 


assuming  the  mean  square  relative  velocities  are  same  in  an  arbitrary  direction  which  implies  yvjx}  —  —  (w/,z)  • 

Saffman  and  Turner  [21]  have  shown  that  only  the  small  error  is  introduced  in  the  collision  frequency  due  to  this 
anomaly.  In  isotropic  turbulence,  the  mean  square  relative  velocity  for  the  viscous  subrange  is  given  by: 

/r.  r.  \  r  (3“3)2  f  i  i  ] 
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(2.44) 


The  mean  square  relative  velocity  for  the  inertial  subrange  is  given  by: 
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Accordingly,  the  mean  square  relative  velocity  (w-w)  can  be  obtained  using  the  equations  given  above  for  (wi-wi 


and  (ws  •  nv ) .  For  inertial  subrange,  the  term  (ws  -ws)  is  neglected  as  explained  in  the  next  section. 


2. 3. 1.4  Collision  frequency  function 

The  collisions  are  likely  to  be  dependent  on  the  relative  motion  between  two  particles.  The  mechanisms  considered 
for  the  collisions  are  the  shear  mechanism  and  the  accelerating  (inertia)  mechanism.  The  shear  mechanism  is  due  to 
the  relative  motion  induced  by  the  viscous  force  inside  the  turbulent  eddy.  The  accelerating  mechanism  is  due  to  the 
relative  motion  induced  by  inertial  effects  between  the  drop  and  suspending  fluid.  The  most  widely  used  collision 
frequency  model  considering  the  shear  mechanism  only  is  the  Saffman  and  Turner  [2 l]’s  model  in  viscous  subrange 
which  is  as  follows: 

P„  =1.294  -  Rl  (2.46) 

\v) 

Saffman  and  Turner  [21]  have  also  derived  the  collision  frequency  including  the  accelerating  mechanism.  However,  as 
explained  in  the  previous  section,  it  is  assumed  that  the  time  scale  of  each  drop  is  smaller  than  the  Kolmogorov’s  time 
scale  (time  scale  of  the  smallest  eddy)  fundamentally  and  this  is  not  matched  with  the  aluminum  drops  in  the 
combustion  chamber  due  to  pp  »  pg.  Instead,  the  collision  kernel  is  modified  to  include  Williams  and  Crane  [24] ’s 
result  for  small  drops  which  is  given  in  the  previous  chapter. 

For  isotropic  turbulence,  the  collision  kernel  which  is  obtained  in  the  previous  chapter  for  the  spherical 
formulation  is  given  by 

f3t=2^R2c{\wr\)  (2.47) 

Using  the  inertial  and  velocity  gradient  terms  given  in  the  previous  chapter,  the  mean  square  relative  velocity  for 
viscous  subrange  is  given  as  follows: 
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The  mean  value  of  radial  relative  velocity  (|wr|^  is  independent  of  the  orientation  of  radial  direction  in  isotropic 
turbulence.  In  analogy  with  Saffman  and  Turner  [21],  it  is  assumed  that  wr  is  aligned  with  the  x-axis  so  that 
^ | wr  | —  (|vt’t|)  .  Assuming  the  mean  square  relative  velocities  are  same  in  an  arbitrary  direction  which  implies 


wr 


~{wy) ~iw^) '  Again,  this  anomaly  leads  to  only  small  error  as  discussed  in  the  previous  chapter.  Thus,  the 


mean  square  relative  velocity  in  radial  direction  is  given  by 
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In  analogy  with  Williams  and  Crane  [24],  it  is  assumed  that  the  relative  velocity  in  radial  direction  spreads  by  a 
Gaussian  distribution.  Thus,  the  mean  of  absolute  relative  velocity  is  the  first  order  moments  of  |u';.  which  is  given 
by 
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Accordingly,  the  collision  frequency  function  for  the  viscous  subrange  is  given  by 
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The  same  discrepancy  as  Saffman  and  Turner  [21]’s  two  models  is  observed.  When  the  drops  are  identical  (inertial 
effects  of  two  adjacent  drops  confined  in  the  smallest  eddy  are  same),  the  constant  becomes  1.671  whereas  the 
constant  in  the  model  concerned  the  shear  mechanism  only  is  1.294.  This  discrepancy  is  caused  by  the  different 
approximation  in  defining  isotropy  as  explained  by  Saffman  and  Turner  and  the  error  is  considered  as  small  [21]. 

For  inertial  subrange,  Kuris  and  Kusters  [31]  have  shown  that  the  accelerative  mechanism  becomes  more 
dominant  with  increasing  particle  sizes  comparing  Saffman  and  Turner’s  shear  term  and  the  universal  solution  of 
inertial  term  derived  in  their  analysis  which  is  very  similar  to  Williams  and  Crane’s  solution.  The  mean  square  of  the 
velocity  gradient  in  inertial  subrange  is  given  by  [32] 
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Comparing  the  shear  term  in  viscous  and  inertial  subrange,  Saffman  and  Turner  [21] ’s  term  is  second  order  in  particle 
size  whereas  Rotta  [32] ’s  term  is  <9(i?2/3) .  Therefore,  the  effect  of  velocity  gradient  in  inertial  subrange  is  neglected 
and  the  mean  square  relative  velocity  in  radial  direction  and  the  collision  frequency  function  for  the  inertial  subrange 
are  given  as  follows: 
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(2.53a) 


(2.53b) 


2. 3. 1.5  Collision  frequency  kernel  -  Laminar  flow 

The  first  consideration  of  particle  collision  rate  in  a  laminar  flow  goes  back  to  1917  by  Smoluchowski  [13].  He 
considered  a  uniform  laminar  shear  flow,  which  is  a  special  case  in  which  only  one  directional  velocity  component  of 
the  fluid  exists  and  it  varies  linearly.  He  assumed  that  the  particles  follow  the  exact  flow  streamlines  so  that  the  shear 


rate,  G,  of  the  flow  motion  is  responsible  for  the  relative  motion  between  two  particles.  He  also  assumed  that  the 
streamlines  are  straight  so  that  the  effect  of  the  particle  on  the  flow  motion  is  neglected  in  his  analysis.  Considering 
two  particles  in  a  laminar  flow  which  are  a  fixed  central  particle  of  radius  r7  and  flowing  particle  of  radius  rh  the 
collision  occurs  when  a  center  of  flowing  particle  approaches  within  a  distance  (rt  +  rj).  A  distance  (rf  +  ry)  is 
conventionally  defined  as  a  radius  of  collision  sphere.  This  collision  is  illustrated  in  Fig.  2.4a.  The  collision  frequency 
can  then  be  calculated  by  the  flux  of  the  fluid  inwardly  normal  to  the  effective  cross-section  collision  area  (= 

)  multiplied  by  the  number  density  of  moving  particles.  The  corresponding  cylindrical  formulation  of  collision 
frequency  is  described  Fig.  2.4b.  The  radius  of  cylinder  is  the  collision  radius,  (rz  +  r7)  and  the  axis  of  cylinder  passes 
through  the  center  of  the  fixed  particle.  The  total  flux  towards  the  effective  collision  area  (which  is  same  as  the 
collision  frequency  function)  is  given  by 

[il=\>G(ri+rjj  (2.54) 

The  Smoulchowski’s  analysis  is  oversimplified  and  it  is  desirable  to  take  into  account  more  complex  motion  of 
flow  and  accurate  amount  of  flux  across  the  surface  of  a  sphere  in  spherical  formulation  of  collision  frequency 
function.  More  recently,  this  is  first  considered  by  Kramer  and  Clark  [33].  They  obtained  the  orthokinetic  coagulation 
frequency  (or  collision  frequency)  for  the  laminar  incompressible  flow  considering  the  strain  rates  acting  within  the 
fluid.  Although  their  approach  can  be  considered  as  more  comprehensive  than  Smoulchowski’s  analysis,  it  is  still 
limited  because  the  strain  rates  which  cause  departure  of  a  moving  particle  from  a  fixed  central  particle  are  simply 
eliminated  from  the  collision  process.  More  precisely,  the  positive  strain  rate  can  contribute  to  the  amount  of  flux 
flowing  towards  the  surface  of  the  collision  sphere  because  the  negative  and  positive  strain  rate  can  be  applied  at  the 
same  time.  Thus,  the  positive  strain  rate  will  result  in  decrease  of  the  total  flux  and  it  should  not  be  ignored.  In 
addition,  more  practically,  the  compressible  effect  can  deform  the  fluid  element  and  it  can  contribute  to  the  particle 
collision  because  the  contraction  of  element  occurs  ( dUi  /  dx.  <  0 )  and  then  more  contribution  to  the  collision  will 
take  place  than  incompressible  case  ( dUi  /  dx.  =0).  This  effect  can  be  significant  in  case  of  the  flow  in  a  rocket 
chamber  or  nozzle.  Therefore,  their  approach  will  be  modified  here  to  incorporate  the  neglected  effects  mentioned 
above. 


A? 


Fig.  2.4  Schematic  of  Smoulchowski’s  droplet  collision;  (a)  illustrates  two  unequal 
droplets  in  uniform  laminar  flow,  (b)  shows  the  cylindrical  formulation  of 
collision  process 


Considering  the  motion  of  fluid  elements  of  viscous  fluid,  it  can  be  thought  that  the  fluid  element  responds  linearly  to 
deformation  rate  (or  velocity  gradient,  deformation  rate  tensor  is  etj  =  dUi  /  dxj )  when  the  fluid  element  is  small.  In 


addition,  the  values  of  deformation  rate  acting  on  the  fluid  element  can  be  considered  as  constants  as  long  as  the  fluid 
element  is  sufficiently  small.  When  two  adjacent  particles  are  assumed  to  exist  in  the  fluid  element,  the  motion  of  the 
particles  is  determined  by  the  deformation  of  fluid  element.  When  the  fluid  element  deforms  due  to  the  velocity 
gradient,  the  fluid  elements  experience  the  linear  deformation  which  can  cause  the  volume  change  of  element  in  a 
compressible  fluid,  the  rotation,  and  angular  deformation  which  change  the  shape  of  the  fluid  element.  The  types  of 
fluid  element  deformation  are  illustrated  in  Fig  2.5. 

The  deformation  rate  tensor,  eih  can  be  represented  by  a  linear  combination  of  two  2nd  rank  tensors  as  follows: 
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The  first  term  in  RHS  is  antisymmetric  part  of  the  deformation  rate  tensor  which  is  called  as  the  rate  of  rotation  tensor, 
Q ip  and  the  second  term  is  symmetric  part  called  as  the  rate  of  strain  tensor,  sip  which  are  given  by: 
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(2.56b) 


The  rate  of  rotation  tensor  is  responsible  to  the  rotation  of  the  fluid  element.  Considering  the  two  particles  in  the 
fluid  element,  the  rotation  does  not  induce  the  collision  as  described  in  Fig.  2.6a.  The  rate  of  rotation  tensor  causes  no 
change  of  shape  or  volume  of  fluid  element  and  correspondingly  it  does  not  cause  the  particle  motion  in  the  direction 
of  reduction  of  the  distance.  The  Fig.2.6b  shows  the  collision  occurred  by  the  normal  components  of  the  rate  of  strain 
tensor  and  Fig.  2.6c  shows  the  collision  occurred  by  the  shearing  components  of  the  tensor.  Therefore,  the  rate  of 
strain  tensor  can  be  thought  as  the  only  source  which  is  responsible  for  the  collision. 
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Fig.  2.5  Schematic  of  the  fluid  element  deformation 
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Fig.  2.6  Schematic  of  influence  on  the  collision  due  to;  (a)  rotation,  (b)  linear 
deformation,  and  (c)  angular  deformation 


Fig.  2.7  Schematic  of  relative  velocity  components  induced 
by  normal  strain  rate 


Following  Kramer  and  Clark  [33]  and  Clark  [34],  the  rate  of  strain  tensor  can  be  diagonalized  without  loss  of 
information  by  the  rotation  of  the  coordinate  system  to  principal  coordinate  because  the  rate  of  strain  tensor  is 
symmetrical.  In  case  of  two-dimensional  symmetric  rate  of  strain  tensor  (^12=^21)  can  be  diagonalized  as  follows: 
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where 


sj  j  =  sl  j  cos2  9  +  s22  sin2  9  +  2 sl2  sin  9  cos  9 
s n  =  sx  i  sin2  9  +  s22  cos2  6  —  2 sl2  sin  9  cos  9 
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Now,  the  collision  is  induced  only  by  the  normal  components  of  transformed  rate  of  strain  tensor  and  the  next  step 
will  be  to  calculate  the  total  flux  across  the  surface  of  the  collision  sphere.  The  Fig.  2.7  illustrates  the  components  of 
relative  velocities  induced  by  normal  strain  rates  in  two  dimensions  with  a  negative  component  in  principal  x 
direction  and  a  positive  component  in  principal  y  direction.  The  velocity  components  induced  by  normal  strains  are 
given  by 

Wx  =  snR  sin  (j) 

(2.58) 

Wy  =*22R,  cos  (f> 

Therefore,  the  velocity  induced  by  normal  strains  in  radial  direction  is  given  by 

Wr  =  s\  {RC  sin 2  ^  +  s'22Rc  cos2  (/)  (2.59) 

Based  on  the  equations  given  above,  the  hydrodynamic  collision  frequency  function  in  laminar  flow  that  is  same  as  the 
inwardly  normal  flux  across  the  collision  sphere  can  be  calculated  as  follows: 

A,,  =•/,/,=-  j  WrdA 
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The  aerodynamic  collision  frequency  function  in  laminar  flow  field  can  be  obtained  from  the  fast  Eulerian 
approach.  From  the  fast  Eulerian  method,  the  relative  velocity  between  two  particles  due  to  slip  motion  assuming  the 

same  DU  /  Dt  at  two  close  points  is  given  by 

DU 


Wr  =(ti—T2)- 


Dt 


(2.61) 


Based  on  the  conservation  of  flux  across  the  particle  surface  ( J  W ,  •  n  dA  =  0 ),  the  inwardly  normal  flux  is  given  by 

7,a=-\  WrdA  =  U\Wr\dA 


wr<  0 

where  W.  =  Wr  •  n 

The  coordinate  transformation  gives  the  following  aerodynamic  collision  frequency  in  laminar  flows: 


(2.62) 


(2.63) 


2.3.2  Breakup  frequency  kernel 

2.3.2. 1  Turbulent  flow  -  Hydrodynamic  breakup 

The  viscous  shear  forces  in  a  turbulent  suspension  acts  on  the  droplet  surface  and  results  in  the  velocity  gradient. 
This  velocity  gradient  leads  to  deform  the  droplet  surface  and  the  breakup  of  the  drop  may  occur  [35,  36].  Considering 
the  hydrodynamic  stresses  as  a  source  of  the  breakup,  Delichatsios  and  Probstein  [37]  have  derived  the  breakup 
frequency  in  the  inertial  subrange.  They  have  used  the  approximation  of  the  probability  density  distribution  of  that  the 

velocity  difference,  A  w(A),  across  the  drop  of  diameter  D,  to  Gaussian  distribution  with  the  variance  ^(A  u)2^  and  cut¬ 


off  velocity  A uc.  The  breakup  of  a  droplet  occurs  when  the  velocity  difference  exceeds  its  critical  value  A ucrit.  The 
breakup  frequency  function  (breakup  frequency  kernel)  b(R of  the  droplet  of  radius  Rf  and  the  distribution  function 
are  given  by 
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respectively.  In  their  model,  the  variance  cr2  and  the  cut-off  velocity  is  given  by 

V 

(72=((Am)2)  =  1.88(2^)2/3 
Amc=3((Am)2)  =3ctv 


(2.65) 


respectively.  Regardless  of  the  direction  of  velocity  difference  acting  on  a  droplet,  all  of  the  velocity  differences  can 
be  responsible  for  the  breakup.  Thus,  the  estimation  of  the  velocity  difference  A u  is  done  in  the  way  of  finding  the  first 
order  moment  of  |Aw|  instead  of  using  the  root-mean-square  of  the  velocity  difference.  However,  the  first  order 

moment  of  the  above  probability  distribution  on  |Aw|  does  not  converge.  Instead,  we  assume  that  the  velocity 
difference  spreads  by  the  exact  Gaussian  distribution  (no  term  of  the  cut-off  velocity,  in  analogy  with  Williams  and 
Crane  [24]).  Thus,  the  velocity  difference  assuming  it  is  the  first  order  moments  of  |Aw|  is  given  by 


Au  =  J  \Au\P(Au)d(Au) 
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where  P(Au)  is  the  Gaussian  distribution.  Therefore,  applying  the  above  equation  into  the  given  distribution  function, 
the  breakup  frequency  for  inertial  subrange  is  expressed  by 
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It  should  be  note  that  the  probability  distribution  P(Au )  is  negative  when  8crit/8  is  smaller  than  53.154  due  to  the 
existence  of  the  term  of  the  cut-off  velocity.  Thus,  the  breakup  frequency  is  set  to  zero  in  this  limit. 

Kusters  [38]  have  derived  the  breakup  frequency  by  assuming  that  the  velocity  difference  A u(Dt)  across  the  drop 
of  diameter  Dt  follows  the  exact  Gaussian  distribution  for  the  viscous  and  inertial  subrange  (in  analogy  with  Saffman 
and  Turner  [21]  and  Williams  and  Crane  [24]).  His  work  has  started  from  the  same  formula  of  the  breakup  frequency 
given  by  Delichatsios  and  Probstein  [37]  except  the  assumption  of  the  probability  distribution.  We  used  his  approach 
in  modeling  of  the  breakup  frequency  in  viscous  subrange.  The  breakup  frequency  and  Gaussian  distribution  are  given 

by 
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The  variance  a1  and  the  mean  velocity  difference  assuming  it  is  the  first  order  moments  of  |Am|  are  given  by 
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Therefore,  the  breakup  frequency  for  viscous  subrange  is  given  by 
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The  hydrodynamic  stress  on  a  droplet  surface  due  to  the  viscous  force  can  be  characterized  by  Capillary  number, 
Ca,  which  is  the  ratio  of  the  stress  due  to  the  continuous  phase  velocity  gradient  to  the  stress  on  the  droplet  surface. 
The  Capillary  number  is  usually  used  in  the  analysis  of  the  viscous  force  acting  on  the  drop  without  the  inertial  force. 
The  classical  definition  of  the  capillary  number  on  a  droplet  radius  Rf  is  given  by 

mgr, 


Ca  =  - 


(2.71) 


where  //  is  the  suspension  flow  viscosity,  G  is  the  velocity  gradient,  and  a  is  the  surface  tension.  The  breakup  can  be 
thought  as  it  is  occurred  when  the  Capillary  number  exceeds  the  critical  Capillary  number  Cacrit  [39].  The  critical 
Capillary  number  depends  on  the  flow  type  and  the  viscosity  ratio  between  a  drop  and  suspension  flow  fip/  ju  [40].  The 
velocity  gradients  induced  only  by  the  viscous  force  in  viscous  and  inertial  subrange  is  given  by 
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respectively.  Therefore,  the  critical  energy  dissipation  rate  corresponding  to  the  critical  Capillary  number  in  viscous 
and  inertial  subrange  is  given  by 
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respectively.  Due  to  the  lack  of  experimental  data,  the  critical  Capillary  number  is  obtained  from  the  results  of 
immiscible  fluids  experiment  in  a  simple  shear  flow.  Cristini  et  al.  [41]’s  numerical  result  gives  that  the  critical 
Capillary  number  is  approximately  0.46  at  the  1.3  viscosity  ratio  which  approximates  the  condition  inside  the 
combustion  chamber. 


2. 3. 2. 2  Turbulent  flow  -  Aerodynamic  breakup 

The  breakup  mechanism  of  liquid  drops  in  a  gas  suspension  is  usually  characterized  by  the  aerodynamic  forces  (or 
inertial  forces)  based  on  the  relative  velocity  between  the  gas  and  droplet.  The  non-dimensional  parameters  used  in  the 
breakup  due  to  aerodynamic  forces  are  the  Weber  number  and  Ohnesorge  number  given  as  follows: 

i2 

A  Ll 

! -  Oh=  .  p  (2.74) 

The  Weber  number  is  the  ratio  of  the  fluid  inertia  to  the  surface  tension  and  Ohnesorge  number  is  the  relative 
importance  of  the  viscous  forces  to  inertial  and  surface  tension  forces.  Here,  | U  —  UP,^  is  the  velocity  difference 

between  the  droplet  i  and  surrounding.  The  degree  of  deformation  and  breakup  is  determined  by  these  parameters. 
According  to  Hsiang  and  Faeth  [42],  there  is  no  breakup  observed  when  the  Ohnesorge  number  is  larger  than  4. 
Because  this  is  not  our  case  (the  Ohnesorge  number  of  A1/A1203  particles  in  a  chamber  is  typically  smaller  than  4 
under  high  temperature  condition),  the  Weber  number  becomes  the  only  parameter  relating  with  breakup.  The  breakup 
occurs  when  the  droplet  Weber  number  is  larger  than  the  critical  Weber  number.  Thus,  the  slip  velocity, 

A u  =  I U  —  Up,i  I ,  corresponding  to  the  critical  Weber  number  is  given  by 
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The  aerodynamic  breakup  frequency  for  both  of  viscous  and  inertial  subrange  is  assumed  to  be  given  by  the 
following  using  the  aerodynamic  particle  break  time  and  the  breakup  efficiency  corresponding  to  the  critical  velocity 
difference: 
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where  the  constant  C  is  obtained  from  the  empirical  correlation  for  the  breakup  time  given  by  Hsiang  and  Faeth  [42], 
The  constant  C  is  given  by 
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for  We  <  103,  0.1  <Oh  <3.5 


Following  Levins  and  Glastonbury  [28],  the  slip  velocity  A us  can  be  related  to  turbulence  parameter.  They  have 

started  from  Tchen  [27] ’s  force  balance  equation  between  a  drop  and  surroundings  and  the  applicable  drops  are  not 
limited  by  their  relative  time  scale  to  the  time  scale  of  the  smallest  eddies.  In  contrast  to  Williams  and  Crane  [24]  who 
used  a  Stokesian  particle  assumption,  Levins  and  Glastonbury  [28]  used  actual  drag  coefficients  corresponding  to  drop 

size  classes  can  be  used  here.  Following  their  approach,  the  mean  square  slip  velocity  q1  =  ({u  -  U p  j  •  (u  -  U p  for  a 

random  turbulent  fluctuation  case  is  given  as  follows  assuming  the  exponential  form  of  the  Lagrangian  correlation 
function  [28]: 
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Here,  L  is  the  integral  length  scale  of  the  largest  energy-containing  eddy  defined  in  above  chapter.  The  root  mean 


square  fluctuating  velocity  (u2^ 


1/2 


is  given  by  (2k  /  3)12  in  isotropic  turbulence.  Therefore,  the  mean  square  slip 


velocity  q  can  be  obtained  by  solving  the  above  equation  numerically.  Assuming  the  Gaussian  distribution  of  relative 
velocity,  the  slip  velocity  A u  can  be  given  as  (21  n\12  q  by  calculating  the  first  order  moment  of  C/-C/J  . 
Newton’s  method  or  bisection  method  is  used  here.  The  drag  coefficient  is  given  by  [43] 
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where  Re  =  pAu  Dip.  The  appropriate  critical  Weber  number  for  A1203  is  found  in  Caveny  and  Gany  [43]  which 

p  si 

is  given  as  28.  The  larger  breakup  frequency  is  chosen  comparing  the  breakup  frequencies  induced  by  viscous  and 
inertial  force  (comparing  hydrodynamic  and  aerodynamic  breakup)  for  both  of  viscous  and  inertial  turbulence 
subranges. 


2. 3. 2. 3  Laminar  flow  -  Hydrodynamic  breakup 

Following  the  same  approach  given  in  hydrodynamic  collision  frequency  model  in  laminar  flow,  the  average 
velocity  difference  across  a  droplet  of  radius  using  the  Euclidean  norm  of  rate  of  stress  tensor  is  given  by. 

AC/  =  Ri^S  sf]  +  24  +  s2}  (2.80) 

Once  the  velocity  difference  AU  across  the  droplet  is  obtained,  this  can  be  used  in  the  hydrodynamic  breakup 
model.  The  hydrodynamic  breakup  model  in  laminar  flow  is  given  by 
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The  critical  velocity  difference  across  the  droplet  A  Ucrit  can  be  obtained  in  terms  of  the  critical  Capillary  number 
Cacrit  as  follows: 
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(2.82) 
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2. 3. 2.4  Laminar  flow  -  Aerodynamic  breakup 

Using  a  fast  (or  equilibrium)  Eulerian  approach,  the  velocity  of  particle  phase  can  be  obtained  then  the  mass  and 
momentum  conservation  are  no  longer  needed  to  be  calculated.  Fundamentally,  this  will  allow  us  to  consider  slip 
between  phases  in  the  model  and  the  resultant  aerodynamic  forces  on  droplets  that  may  or  may  not  lead  to  breakup.  In 
the  fast  Eulerian  approach,  the  particle  phase  velocity  is  handled  as  a  field  variable  with  a  limitation  of  small  size  and 
large  density  of  particles  and  the  relative  velocity  between  the  particle  and  the  surrounding  fluid  is  given  by 
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where  U p  and  U  are  the  particle  phase  and  gaseous  phase  velocity  vectors,  respectively.  The  term  r  is  the  relaxation 
time  of  the  particle  and  D  /  Dt  is  the  material  derivative  in  Eulerian  view  point.  In  analogy  with  Saffman  and  Turner 
[21],  assuming  that  the  carrier  fluid  velocities  near  two  adjacent  particles  are  same,  the  velocity  difference  of  the 
carrier  fluid  between  two  close  point  is  neglected  (this  velocity  difference  is  already  considered  in  the  previous  chapter) 
and  then  the  relative  slip  velocity  between  two  particles  is  given  by 
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Once  the  slip  velocity  ( A U  =  |g| )  is  obtained,  the  modeling  of  aerodynamic  breakup  in  laminar  flow  can  be 


considered  using  the  aerodynamic  breakup  frequency  function: 
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where  the  constant  C  and  critical  slip  velocity  are  given  by  [42] 
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in  terms  of  critical  Weber  number,  Wecrit. 


2.3.3  Coalescence  efficiency 

As  Salita  [44]  explained,  the  assumption  that  all  collisions  result  in  coalescence  has  usually  been  used  in  the 
multiphase  simulation  of  rocket  combustion  chamber.  However,  many  prior  studies  (i.e.,  Ashgriz  and  Poo  [45])  have 
shown  experimentally  that  all  collisions  of  two  water  drops  does  not  result  in  coalescence.  Four  different  types  of 
collision  outcomes  were  observed  in  the  experiments  as  it  is  explained  previously:  bouncing,  permanent  coalescence, 
temporary  coalescence  followed  by  separation  and  temporary  coalescence  followed  by  a  set  of  satellite  drops.  In 
bouncing  mode,  the  drops  are  reflected  in  the  reverse  direction  of  their  approaching  velocity  due  to  the  effect  of  the 
fluid  film  between  the  drops.  The  permanent  coalescence  refers  to  the  creation  of  stable  drops  after  the  drops  are 
merged.  The  temporary  coalescence  with  separation  is  when  the  merged  drop  is  unstable  and  the  merged  drop  is 
separated  into  two  or  more  drops.  The  temporary  coalescence  with  satellite  drops  is  similar  to  separation  mode  but  it  is 
disintegrated  into  a  set  of  very  small  satellite  drops.  Although  the  last  three  modes  concern  the  phenomenon  after  the 
drops  are  contacted  each  other,  the  first  bouncing  mode  concerns  the  drops  before  they  are  contacted. 

Due  to  the  absence  of  the  coalescence  model  of  metal  droplets  (A1  or  A1203  in  the  rocket  propulsion),  Salita  [44] 
was  motivated  to  perform  a  series  of  coalescence  experiments  using  mercury  drops,  whose  the  density  is  13.5  times 
and  surface  tension  6.5  times  bigger  than  water  drops  at  room  temperature.  By  using  the  mercury  drops,  they  can 
provide  surface  tension  values  near  that  of  A1203.  To  compare  the  coalescence  model  with  the  experimental  results, 
they  have  used  the  water  drops  coalescence  model  of  Brazier-Smith  et  al.  [46].  The  model  of  Brazier-Smith  et  al.  [46] 
postulates  that  the  collision  of  drops  always  result  in  the  unstable  coalescence  and  then  the  merged  drop  will  be 
separated  into  the  same  size  of  incident  drops  if  the  rotational  energy  of  the  merged  drop  exceeds  the  surface  energy 
holding  it  together  [44].  They  concluded  that  the  coalescence  model  of  water  drops  accurately  predicted  that  of 
mercury  drops.  The  coalescence  model  of  water  drops  is  used  here. 

The  processes  of  permanent  coalescence  and  disintegration  are  described  in  Figure  2.8.  Considering  two  particles 
of  radii  rf  and  r7,  the  moving  particle  is  the  particle  of  radius  rt  supposing  the  particle  of  radius  r7  as  a  fixed  particle.  A 
temporally  formed  agglomerate  sphere  due  to  the  collision  has  a  mass  m^nij  and  the  corresponding  radius 
Ro  =  (^3  +  Vj  )1/3 .  The  resulting  rotational  energy  of  the  temporal  agglomerate  from  the  induced  angular  momentum  by 


the  impact  of  moving  particle  to  the  fixed  particle  is  given  by 


Er=^Pp 


U  p,rel 


x  r  rj 


R: 


(2.87) 


Here  the  separation  is  assumed  to  occur  in  the  way  that  only  two  product  droplets  which  have  the  same  size  as  the 
original  two  droplets  because  Brazier-Smith  et  al.  [46]’ s  experimental  results  for  water  droplets  have  shown  a  good 
agreement  with  the  modeling  using  this  separation  model.  The  energy  required  to  separate  the  temporal  agglomerate 
into  droplets  of  radii  r  is  given  by 


(2.88) 


Es  =  47 TCJ 


Er*2-Ro 


n—l 


Thus,  in  case  of  separation  into  the  two  droplets  which  have  the  same  size  as  the  original  two  droplets  is  given  by 

Es,s  —  47RJ  +  Vj  —  Rl  (2.89) 

Therefore,  the  permanent  coalescence  occurs  when  ER  <  Es,s,  and  this  condition  gives  that 
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Defining  the  coalescence  efficiency  as  the  ratio  of  collision  cross  section  ( nx2)  to  the  maximum  available  collision 
cross  section  nR 2 ,  the  coalescence  efficiency  is  given  by: 


(  \2 
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Rc 

12 


a 


P+P-Rl  R" 


(2.91) 
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where  pp  is  the  density  of  the  droplet,  and  a  is  the  surface  tension  of  droplet. 
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Fig. 2. 8  Schematic  of  coalescence  and  separation  processes  of  collision  between  two  droplets  of 
radii  rt  and  r. 


2.3.4  Collision  efficiency  (Bouncing  efficiency) 

Following  Coulaloglou  and  Tavlarides  [7]  and  Tsouris  and  Tavlarides  [8],  the  drop  collision  efficiency  can  be 
characterized  by  two  terms;  contact  time  and  collision  time.  The  contact  time  is  the  time  of  two  drops  staying  together 
after  they  collided.  The  collision  time  is  the  time  required  for  the  drainage  of  liquid  or  gas  films  between  two  drops.  If 
two  drops  are  staying  together  after  they  collided  for  enough  time  for  the  film  drainage,  the  collision  of  two  drops 
occurs.  Therefore,  the  collision  occurs  when  the  contact  time  is  larger  than  the  coalescence  time.  The  function  of  the 
coalescence  efficiency  derived  by  Ross  [47]  is  given  by  Tsouris  and  Tavlarides  [8]  as  follows: 
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a  =  —  exp 
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(2.92) 


where  t  is  the  average  contact  time,  T  is  the  average  coalescence  time,  and  o  is  the  standard  deviation  for  the 


coalescence  time.  This  equation  can  be  simplified  by  assuming  that  a  is  zero  as  follows  [8]: 


a 


-  exp 


( 


(2.93) 


Considering  the  liquid  droplets  in  gaseous  fluid  system,  the  deformation  of  the  interface  between  two  droplets  is 
supposed  to  be  significant  unlikely  to  the  solid  particles  in  fluid.  Thus,  assuming  the  deformation  of  the  interface  is 
significant,  the  parallel-film  approach  can  be  used  in  this  system.  The  Figure  2.9  illustrates  the  idealized  deformation 
of  the  interface  in  the  parallel-film  approach. 


Fig. 2. 9  Schematic  of  idealized  collision  process  between  two  deformable 
droplets  of  equal  radius  Req 


[48] 


Based  on  the  classical  lubrication  theory,  the  film  drainage  process  can  be  described  by  two  equations  given  by 


—  ( phr )  H - ( U  phr )  =  0  (continuity  equation) 

dt  dr^  m  ’ 

SP  u 

—  +  12//  — | L  =  0  (r-momentum  equation) 


(2.94) 


dr  '  h 2 

where  p  is  the  gas  density,  p  is  the  dynamic  viscosity,  Um  is  the  gap-averaged  radial  velocity  of  draining  flow,  and  P  is 
the  gap-averaged  pressure.  Considering  incompressible  gas  flow  (p  is  constant),  two  governing  equations  can  be 
solved  for  the  pressure  in  the  gap  between  two  spherical  drops  as  follows: 

P(r,t)  =  -lpr2A  (2.95) 

n 

Neglecting  Van  der  Waals  attraction  and  electrostatic  double -layer  repulsive  force,  the  balancing  of  inertial  force  by 
suspension  flow  and  dynamic  force  exerted  by  the  pressure  inside  the  gap  closes  the  system.  The  force  balance 
equation  is  given  by 


3  wr  h 


-  0 


(2.96) 


2  h2 

where  F  is  the  compressive  inertial  force.  Following  Chesters  [23],  the  amount  of  deformation  is  related  to  the 
compressive  inertial  force  as  follows: 


F  =  nr1 
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(2.97) 


Thus,  eliminating  the  radius  of  the  interfacial  circle  rs  from  the  force  balance  equation,  the  rate  of  film  thickness 
change  is  expressed  by 


dh  8  no1 


dt  3  pR  F 


- h 3 


(2.98) 


Consequently,  assuming  the  constant  force  F,  the  calculation  of  the  time  required  for  film  drainage  when  the  drops 
deform  is  given  by 
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where  the  /?,  and  h0  is  the  final  gap  distance  and  initial  gap  distance  between  particles  respectively.  This  collision  time 
is  also  used  in  Coulaloglou  and  Tavlaride  [7], 

Here,  F  is  the  force  acting  on  either  ends  of  two  adjacent  drops.  By  assuming  the  force  is  propotional  to  the  mean 
square  velocity  difference  at  either  ends  of  the  drops,  Coulaloglou  and  Tavlarides  [7]  have  used  the  following 
expression  for  turbulent  flow: 


F  =  pR 2  /  Am  2  (27?  )\ 

eq  \  c  / 


(2.100) 


where  ^ Au2(2R  )^j  is  the  velocity  difference  across  the  distance  2 Rc  in  the  turbulent  flow.  They  have  assumed  that  the 

proportionality  constant  in  this  equation  is  unity  following  Rotta  [32].  Similarly,  the  following  equation  for  the  force 
acting  on  either  ends  of  two  adjacent  drops  (at  a  distance  of  2Rc)  in  laminar  flow  is  used: 

F  =  pR2AU2(2R  )  (2.101) 

eq  c 

where  the  average  velocity  difference  is  obtained  from  Eq.  2.79  at  a  distance  of  2 Rc.  The  velocity  differences  in 
viscous  subrange  and  inertial  subrange  are  given  by 


(2.102a) 
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respectively.  Therefore,  F  for  viscous  subrange  and  inertial  subrange  are  given  by 
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(2.103a) 

(2.103b) 


respectively.  The  final  gap  distance  hf  between  particles  which  is  set  as  300  armstrong  meters  which  is  the  minimum 
of  reasonable  film  thickness  appearing  in  Coulaloglou  and  Tavlarides  [7]  and  the  initial  gap  distance  h0  is  set  as  0.17?^ 
following  Tsouris  and  Tavlarides  [8].  The  contact  time  for  the  viscous  subrange  and  inertial  subrange  is  set  as 
Kolmogorov’s  time  scale  and  Taylor’s  micro  time  scale,  respectively.  The  contact  time  for  the  viscous  subrange  and 
inertial  subrange  is  given  by 
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respectively.  Following  Chesters  [23],  the  contact  time  for  laminar  flow  is  given  as  follows: 


t  = 
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(2.105) 


2.3.5  Implementation  of  the  mean  velocity  effect  in  turbulent  flow 

The  mean  velocity  effect  in  turbulent  flow  on  the  collision  process  can  be  obtained  by  using  decomposition  of  the 
relative  velocity  between  two  particles  which  is  suggested  by  Saffman  and  Turner  [21]  as  follows: 

(^2}  =  (e2-Si  2^  +  {u2-U\  2)  (2.106) 

Here,  the  square  means  the  dot  product.  Using  the  Reynold’s  decomposition,  the  mean  square  relative  velocity 
between  particles  becomes 

(^2}=  (e2)-(ei)2+(72-71  2)+  (u2)-(uy  +(u2-u,  2)  (2.107) 

The  lower  cases  q  and  u  are  fluctuating  parts  of  the  slip  velocity  and  the  surrounding  fluid  velocity.  The  first  and 
third  terms  on  the  RHS  are  given  by  Eqs.  2.61  and  2.79  respectively  (the  average  velocity  difference  across  at  a 
distance  of  Rc  is  obtained  from  Eq.  2.79).  The  sum  of  second  and  fourth  terms  is  obtained  from  Eq.  2.49  for  viscous 
subrange  and  Eq,  2.53a  for  inertial  subrange  (the  fourth  term  is  neglected  in  inertial  subrange).  This  velocity 
difference  is  used  to  calculate  the  coalescence  efficiency  in  turbulent  including  significant  mean  flow  effects  by 
putting  the  above  equation  into  Eqs.  2.90.  In  case  of  collision  efficiency,  it  is  assumed  that  the  force  is  propotional  to 
the  mean  square  velocity  difference  at  either  centers  of  the  drops.  In  addition,  the  minimum  value  between  the 
collision  and  coalescence  efficiencies. 

Also,  the  calculation  of  mean  square  velocity  difference  across  a  droplet  ^(A U)2^  and  the  mean  square  slip 

velocity  U  j  ^  is  important.  Relaxing  the  assumption  of  constant  mean  velocity  and  using  Reynold’s 
decomposition,  these  terms  can  be  given  by 


respectively. 
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(2.108) 


Additionally,  the  key  parts  of  implementation  of  mean  velocity  effect  in  collision/breakup  process  are  the 
calculation  of  collision  and  breakup  frequency  function.  It  is  commonly  assumed  that  the  collision  mechanisms  due  to 
Brownian  motion,  fluid  shear,  and  differential  sedimentation  are  independent  from  each  other  and  the  collision 
frequency  functions  are  additive  [49].  The  collision  induced  by  differential  sedimentation  is  very  similar  to 
aerodynamic  collision  because  the  source  of  the  collision  is  the  velocity  difference  between  two  drops  induced  by 
different  response  of  each  drops  to  their  surrounding  fluid.  It  is  further  assumed  here  that  the  collisions  induced  by 
laminar  and  turbulent  flow  are  independent.  Thus,  the  sum  of  collision  frequency  functions  of  laminar  hydrodynamic 
collision,  laminar  aerodynamic  collision,  and  turbulent  collision  is  used  for  the  total  collision  frequency  function.  In 
case  of  breakup  frequency,  the  maximum  value  between  laminar  hydrodynamic  breakup,  laminar  aerodynamic 
breakup,  and  turbulent  breakup  is  used  for  the  calculation.  The  summary  of  equations  used  for  collision  and 
coalescence  efficiencies,  and  collision  and  breakup  frequency  function  is  provided  in  Table  2.1. 


Table  2.1  Summary  of  equations  used  for  collision/coalescence  efficiencies  and  collision/breakup  frequency  functions 


Laminar 

Turbulent 

Turbulent  with 
mean  flow  effects 

Collision 

efficiency 

Eq.  2.93 

-  Eqs.  2.99,  2.101,  and  Eq. 
2.105  are  used 

Eq.  2.92 

Eq.  2.93 
-  the  velocity 
difference  is 
obtained  from  Eq. 
2.108 

-  the  contact  time  is 
assumed  as  the  sum 
of  laminar  and 
turbulent 

Viscous 

Inertial 

Eqs.  2.98,  2.103a, 
and  2. 104a  are  used 

Eqs.  2.99,2.103b, 
and  2. 104b  are  used 

Coalescence 

efficiency 

Eq.  2.91 

-  the  amount  of  relative  velocity 
is  obtained  from  the  sum  of  Eq. 
2.80  at  a  distance  of  Rc  and  the 
absolute  value  of  Eq.  2.61 

Eq.  2.91 

Eq.  2.91 

-the  relative  velocity 
is  obtained  from  Eq. 
2.107 

Viscous 

Inertial 

Eq.  2.49  used 

Eq.  2.53a  used 

Collision 

frequency 

function 

Hydrodynamic 

Aerodynamic 

Viscous 

Inertial 

Sum  of  Eq.  2.60, 
2.63,  and  2.51  or 
2.53b 

Eq.  2.60 

Eq.  2.63 

Eq.  2.51 

Eq.  2.53b 

Breakup 

frequency 

function 

Hydrodynamic 

Aerodynamic 

Hydrodynamic 

Aerodynamic 

Maximum  between 
Eqs.  2.81,  2.85, 
2.70,  2.67,  and  2.76 

Eq.  2.81 

Eq.  2.85 

Viscous 

Inertial 

Eq.  2.76 

Eq.  2.70 
-Eq.  2.73a 
used 

Eq.  2.67 
-Eq.  2.73b 
used 

III.  Results 


3.1  Problem  Statement 

3.1.1  Baseline  geometry 

A  series  of  simulation  was  performed  to  validate  the  code  comparing  the  results  to  Crowe  et  al.[50]’s  experimental 
results  and  Hermsen’s  correlation.  Crowe  et  al.’s  nozzle  has  a  length  of  5.2  cm  and  a  inlet  radius  of  1.587  cm.  The 
throat  is  located  at  2.113  cm  downstream  from  the  inlet  and  its  diameter  is  1.27  cm.  The  corresponding  area  ratio  of 
the  throat  to  the  inlet  is  6.246  and  the  ratio  of  the  throat  to  the  exit  is  also  6.246.  The  figure  3.1  illustrates  their  nozzle 
geometries.  They  also  performed  the  experiments  measuring  the  particle  diameter  at  the  chamber  without  the  nozzle 
attached  by  pressurizing  the  collection  chamber  to  a  certain  level.  The  particles  collected  from  the  motors  containing 
16%  Aluminum  particles  had  the  mass  mean  diameter  of  0.74  pm  and  the  standard  deviation  of  0.456  on  pressure  over 
150psi  with  only  slight  variation  of  values.  Because  we  can  obtain  the  parameters  of  their  experiments  required  to 
perform  simulations  from  Crowe  and  Willoughby [52],  we  decide  to  perform  simulations  under  their  experimental 
conditions. 

A  series  of  simulations  was  also  performed  to  compare  the  results  to  the  correlation  under  the  nozzle  configuration 
used  by  Shegal[53]  for  150psi  chamber  pressure.  His  experiments  performed  with  Polyurethane -type  solid  propellants 
containing  12%  aluminum.  The  motor  dimension  is  5  inch  outer  diameter  by  6  inch  long  with  a  circular  port  of  2  inch 
diameter.  The  conical  convergent  nozzle  is  attached  to  the  motor  and  the  chamber  pressure  was  changed  by  adjusting 
the  throat  diameter.  He  obtained  the  particle  size  information  at  the  nozzle  exit  (or  nozzle  throat)  by  firing  motor  into  a 
tank.  He  reported  the  size  data  from  the  particles  attached  to  tank  wall. 

The  particle  size  information  in  the  motor  of  Shegal  can  be  obtained  from  Fein[54].  While  holding  the  chamber 
pressure  150  and  500  psi  by  pressurizing  the  tank  to  the  desired  level,  the  motors  without  the  nozzle  are  fired  into  the 
tank.  The  measured  MMD  of  particles  were  0.79  and  2.39  pm  for  150  and  500  psi  chamber  pressures,  respectively.  At 
these  chamber  pressures,  the  measured  MMDs  at  the  nozzle  exit  were  1.5  and  3.5  pm,  respectively. 

In  order  to  reduce  the  error  due  to  the  fast  Eulerian  assumption,  we  choose  the  nozzle  geometry  for  Shegal’ s  150 

psi  case.  The  large  particle  diameter  gives  the  large  relaxation  time  (z  —  2pp  +  pg  r2  / 9//  ),  then  it  can  lead  the  large 

amount  of  error  in  the  particle  phase  velocity,  U p  .  According  to  Ferry  and  Balachandar[18],  the  particle  phase 
velocity  can  be  obtained  from  the  fast  Eulerian  approach  within  the  reasonable  error  bound  when  the  relaxation  time  is 
less  than  the  fluid  time  scale  which  is  defined  by  the  inverse  of  the  maximum  of  absolute  compressive  strain).  Thus, 
the  difficulty  specific  to  the  fast  Eulerian  approach  arises,  which  is  that  this  approach  can  produce  negative  values  of 
particle  phase  velocity.  To  overcome  this  difficulty,  clipping  has  been  used.  The  fast  Eulerian  approach  can  be 
replaced  by  solving  momentum  equations  for  the  particle  phase  to  obtain  more  accurate  particle  phase  velocity. 
However,  it  is  out  of  our  concerns  for  the  computational  efficiency. 

Because  Shegal  did  not  provide  the  detailed  geometrical  information,  we  performed  the  isentropic  analysis  and 
obtained  the  nozzle  throat  diameter.  To  account  for  the  particle  size  variation  in  the  diverging  section  and  compare  the 
particle  size  data  at  the  nozzle  exit,  the  conical  diverging  section  with  a  1 8  deg  half  angle.  In  addition,  the  nozzle  exit 
diameter  is  chosen  according  to  the  perfect  expansion  assumption  at  the  sea  level  and  it  gave  the  shock  wave  free 
condition  inside  the  diverging  section. 

Figure  3.2  highlights  a  typical  axisymmetric  unstructured  mesh  used  in  the  computations.  The  nozzle  has  a  length 
of  15.4  cm  and  a  inlet  radius  of  6.35  cm.  The  throat  is  located  at  12.6  cm  downstream  from  the  inlet.  The 
corresponding  area  ratio  of  the  throat  to  the  inlet  is  30.57  and  the  ratio  of  the  throat  to  the  exit  is  2.43.  The  inlet 
geometry  is  horizontally  smoothed  out  to  remove  additional  disturbances  caused  by  sharp  geometry.  A  simulation  is 
performed  in  a  typical  axisymmetric  unstructured  mesh  form.  The  physical  time  consumed  during  a  test  case 
simulation  is  about  four  days  using  4  cpus  a  time  stepping  by  maximum  CFL  number  in  the  message  passing 
environments. 


10  mm 
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Figure.  3.2  Schematic  of  the  test  nozzle  geometry  for  Shegal’s  experiment  (265x80  cell  mesh) 


All  our  current  simulations  are  mnning  on  our  current  HPC  (High-performance  computing)  cluster.  The  hardware 
configuration  of  the  HPC  cluster  is  presented  in  Table  3.1.  The  software  used  for  the  calculation  is  presented  in  Table 
3.2. 

A  non-slip  boundary  condition  is  imposed  at  the  nozzle  wall  and  a  pressure  inlet  and  supersonic  outlet  condition  is 
set  as  boundary  conditions.  The  inlet  kinetic  energy  k  and  specific  dissipation  rate  <d  are  set  as  sufficiently  small 
values.  The  gas  mixture  properties  are  summarized  in  Table  3.3.  These  properties  are  obtained  from  the  properties  for 
solid  propellant  rocket  simulations  of  Lupoglazoff  and  Vuillot  [56]  and  Najjar  et  al  [2].  The  inlet  temperature  is 
obtained  from  Fein[54].  Sutherland’s  law  is  used  for  the  viscosity  rather  than  the  constant  viscosity  assumption  with 
the  reference  temperature  Tref  and  Sutherland’s  constant  Sref  given  in  Table  3.3  and  it  is  given  by 
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Table  3.1  Hardware  configuration  of  current  HPC  cluster 


Hardware 

Description 

Total 

Quantity 

CPU 

Head  node 

2.0  GHz  AMD  Opteron  Quad  Processor 

1 

Computing  nodes 

2.4  GHz  AMD  Opteron  Quad  Processor 

56 

Motherboard 

Head  node 

HP  ProLiant  DL385 

1 

Computing  nodes 

HP  ProLiant  DL145 

56 

Memory 

Head  node 

8  GB 

284  GB 

Computing  nodes 

13  nodes  -  8GB,  43  nodes  -  4  GB 

Storage 

Head  node 

1.2  TB 

5.7  TB 

Computing  nodes 

80  GB 

Switch 

HP  ProCurve  5406zl  (10  GB  Interconnection  support) 

1 

HP  ProCurve  5400zl  (10  GB  Interconnection  support) 

1 

Table  3.2  Software  used  for  calculation 


Software 

Description 

Operating  system 

RedHat  Linux  7.2 

Fortran  Compiler 

PGI  Compiler  7.1-2  by  Potland  Group 
Cluster  Development  kit 

MPI 

MPICH2  1.0.5 

T able.  3.3  Gas  mixture  properties  and  pressure  boundary  conditions 


Quantity 

Value 

MW  (kg/kmol) 

27.76 

Cp  (J/kg  •  K) 

2439.04 

Are/ (kg/m  •  s) 

36.0e-05 

T inlet  =  Tref  (K) 

3279 

V/(K) 

120 

Pin  (N/m2) 

according  to  simulation  cases 

P out  (N/m2) 

101325 

3.1.2  Particle  phase  properties  and  boundary  conditions 

The  density  for  the  particle  phase  is  obtained  from  A1203  density  relationship  given  by  Najjar  et  al.  [2]  as  follows: 

p  =5632-1. 1277  (kg/m3)  (3.2) 

P 

where  T  is  in  deg.  K.  Due  to  the  limited  available  data,  a  surface  tension  and  dynamic  viscosity  of  liquid  A1203  are 
obtained  from  Hatch  [57]  for  molten  Aluminum  instead  of  Aluminum  Oxide.  The  surface  tension  and  dynamic 
viscosity  of  molten  Aluminum  are  given  by 

g  —  0.868  —  0.000152(77  —  T  )  (N-nT1) 


m 

u  =  0.000 1492  exp(l  984.5/7’)  (N-s/m2) 

p 

where  Tm  is  the  melting  temperature  of  Aluminum  and  it  is  given  as  933.47  K. 


(3.3) 


The  mass  mean  diameter  (or  Herdan’s  mean  diameter)  D43  can  be  obtained  from  the  weights  and  abscissas  which 
is  calculated  by 
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(3.4) 


The  purpose  of  the  current  study  is  to  simulate  and  investigate  the  coalescence/breakup  processes  in  the  typical 
converging/diverging  nozzle.  The  coalescence/breakup  processes  are  sensitive  to  the  distribution  of  particles.  The  log¬ 
normal  particle  number  distribution  or  exponential  distribution  can  be  used  following  Najjar  et  al.  [2]  and  Fein  [54], 
respectively.  Najjar  et  al.  have  referred  other  researcher’s  finding  of  lognormal  and  bimodal  size  distribution  of 
droplets  entering  the  chamber  from  the  solid  propellant  surface.  Gany  et  al.  [l]’s  experimental  results  of  the 
distribution  of  the  droplets  leaving  the  propellant  surface  is  close  to  a  lognormal  distribution.  The  model  proposed  by 
Fein  [54]  is  the  exponential  distribution  rather  than  lognormal  distribution.  Fein  compared  his  modeling  with  the 
experimental  data  performed  by  Shegal  and  the  good  agreement  between  the  model  and  experiment  is  obtained. 

Thus,  the  lognormal  distribution  is  given  by 
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where  dn  is  the  number  fraction  of  drops  in  a  given  range  dD ,  a  is  the  standard  deviation,  and  Dm  is  the  mean  drop 

s 

diameter.  The  transformed  coordinate  by  l3=v  (where  v  is  the  volume  of  a  drop)  can  be  expressed  as  l=D  according  to 
DQMOM  approach  and  the  distribution  of  /  is  also  lognormal  distribution.  Therefore,  the  moments  are  given  by 
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The  analytical  expressions  for  the  raw  moments  of  lognormal  distribution  are  given  by 


m  =n  exp(zlnZ)  +  z2cr2/  2) 

z  total  \  m  s  ) 


z  =  0 . 2N- 1 


The  exponential  distribution  modeled  by  Fein  [54]  is  particle  volume  distribution  which  is  given  by 
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Vm  (6v) 

where  Vn  is  the  number  average  particle  volume.  Therefore,  the  volume  distribution  needs  to  be  converted  to  size 
distribution.  Assuming  the  spherical  droplet,  the  size  distribution  is  obtained  as  follows: 

/(/)  -  f(D)  =  3a  D2  f(v) 

r  .  ,  y/3  (3-9) 

=  /lexp [-/LD]  where  A  =  Itt/V  j 

where  A  is  the  rate  parameter  and  av  is  the  shape  factor  for  the  sphere  which  is  given  by  n/6.  The  analytical  expressions 
for  the  raw  moments  of  exponential  distribution  are  given  by 

(z!// tz)  z  =  0 . 2N-1  (3.10) 
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The  raw  moments  of  inlet  particle  distribution  are  used  to  find  weights  ( iy. )  and  abscissas  ( /. )  from  the  PD  algorithm. 

To  use  these  two  distributions  as  an  inlet  condition,  two  and  one  variables  should  be  known  for  lognormal  and 
exponential  distributions,  respectively.  As  described  above,  we  chose  Shegal’s  nozzle  configuration  used  for  150psi 
chamber  pressure.  Thus,  we  used  his  experimental  results  for  particles  in  the  chamber  as  a  nozzle  inlet  condition. 
According  to  Fein[54],  while  holding  the  chamber  pressure  150  psi  by  pressurizing  the  tank,  the  motor  without  the 
nozzle  gives  that  the  rate  parameter  A  is  4858099.849  1/m  and  D43  is  0.79  pm.  However,  due  to  the  slight  disturbance 
occurring  in  PD  algorithm,  the  calculated  value  of  D43  using  this  rate  parameter  is  0.78  pm.  To  remove  this 
undesirable  effect  from  PD  algorithm,  the  rate  parameter  is  slightly  adjusted  and  is  set  as  4804660.751  1/m.  This 
adjustment  will  not  deviate  the  results  significantly.  Finally,  the  particle  concentration  is  given  by  0.24  assuming 
that  all  of  Aluminum  in  the  propellant  is  converted  to  Aluminum  Oxide. 


3.1.3  Experimental  results  by  other  researchers 

Although  the  current  model  has  an  ability  to  predict  the  drop  size  distribution  inside  the  rocket  chamber  and 
nozzle,  it  is  hard  to  validate  the  code  due  to  the  lack  of  experimental  data  of  drop  size  inside  the  chamber.  The  high 
temperature  and  high  velocity  conditions  in  the  rocket  chamber  and  nozzle  make  it  difficult  to  measure  the  particle 
size.  Until  now,  the  experiments  are  performed  to  measure  the  particle  size  at  the  exit  plane  of  nozzle  (i.e. 
Sambamurthi  [58])  and  lots  of  empirical  correlations  are  developed  to  predict  the  particle  size  at  the  exit  as  it  is 
described  by  Hermsen  [10].  Thus,  we  decide  to  validate  the  predicted  particle  size  data  at  the  nozzle  exit  with  the 
empirical  correlation.  The  empirical  correlation  which  is  widely  used  in  the  solid  rocket  industry  is  Hermsen’ s  model 
[10]  for  the  Aluminum  Oxide  particle  size: 

£>43  =  3.6304D0  2932  (l-exp(-0.0008163£Pr  ))  (3.11) 

where  D43  is  the  mass  mean  diameter  (jum),  Dt  is  the  nozzle  throat  diameter  (in.),  4  is  the  particle  concentration  in  the 
chamber  (gmol/100  g),  Pc  is  the  chamber  pressure  (psi),  and  zc  is  the  average  particle  residence  time  in  chamber  (ms). 
The  average  chamber  residence  time  is  given  by 

r  =  p  V  /  m  (3.12) 

c  c  c 

where  pc  is  the  gas  density  in  chamber  (kg/m3),  Vc  is  the  volume  in  chamber  (m3),  and  m  is  the  propellant  mass  flow 
rate  (kg/s).  Because  the  current  DQMOM  modeling  uses  the  total  number  of  particles,  ntotah  instead  of  the  particle 
concentration,  4,  the  total  number  of  particles  can  be  obtained  using  PD  algorithm.  After  the  weights  ( wt )  and 


abscissas  ( )  are  obtained  by  PD  algorithm  for  a  certain  ntotah  a  ntotai  corresponding  to  a  given  4  is  calculate  by  trial 


and  error  method  on  the  following  equation: 
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where  M  is  the  molar  mass  of  A1203  (g/mole)  which  is  given  by  101.94  for  the  current  study.  According  to 

A12°3 

Hermsen[10],  the  standard  deviation  of  the  correlation  is  0.298  and  it  is  corresponding  to  a  deviation  of  D43  of  about 
±35%  due  to  the  data  scatter  obtained  from  various  collection  and  measurement  techniques. 

To  validate  the  modeling,  we  compared  the  predicted  results  to  Crowe  et  al.’s  experiments.  In  addition,  the 
simulation  is  validated  by  comparing  it  with  Hermsen’ s  correlation.  Therfore,  we  constructed  one  test  matrix  for 
Crowe  et  al.’s  experiments  and  three  different  test  matrices  by  varying  the  chamber  pressure,  particle  concentration, 
and,  nozzle  scale  for  Shegal’s  cases.  Tables  3.4-7  summarize  the  inlet  conditions  for  gas  and  particle  phases  and 
nozzle  scales. 


T able.  3.4  The  particle  phase  inlet  boundary  conditions  -  for  Crowe  et  al.  [50] ’s  experiments 


Case  no. 

C-l 

C-2 

C-3  | 

D, 

0.5 

4 

0.277 

Pc 

(chosen  to  compensate  pressure 
variation  according  to  Hermsen  [10]) 

470 

650 

980 

U 

15  1 

W  total 

7.61el5 

10.57el5 

15.93el5  | 

Dm 

0.361 

°s 

0.456 

Table.  3.5  The  particle  phase  inlet  boundary  conditions  -  chamber  pressure  variation 


Table.  3.6  1 

"he  particle  phase  inlet  boundary  conditions  -  particle  concentration  variation 

Case  no. 

M-l 

IY1-2 

IY1-3 

M-4 

INI-5 

M-6  | 

A 

0.904  1 

& 

0.14 

0.16 

0.18 

0.2 

0.22 

0.24  | 

Pc 

550 

4.15 

ft  total 

10.85el5 

12.4el5 

13.95el5 

15.5el5 

17.04el5 

18.6el5  | 

X 

4804660.751  J 

Table.  3.7  The  particle  phase  inlet  boundary  conditions  -  nozzle  scale  variation 


Case  no. 

S-l 

S-2 

S-3 

S-4 

S-5 

S-6 

A 

0.904 

1.808 

3.616 

5.424 

7.232 

9.04 

4 

0.240 

Pc 

550 

4.15 

ft  total 

18.6el5 

X 

4804660.751 

3.2  Grid  convergence  study 

The  grid  convergence  study  was  performed  on  the  geometry  given  in  above  section  for  Shegal’s  150  psi  case.  The 
chamber  pressure  is  150  psi  and  the  total  number  of  particles  at  chamber  is  2.612el5.  The  mass  mean  diameter  at 
chamber  is  0.79  pm  and  the  corresponding  rate  parameter  of  exponential  distribution  is  4804660.751  using  PD 
algorithm.  The  grid  sizes  used  are  132x60  and  265x80.  We  are  also  performing  the  simulation  on  381x100  grids,  but  it 
requires  exceptionally  long  time  to  get  the  converged  solution  (approximately  30  days)  so  it  is  not  included  here.  The 
results  of  grid  convergence  study  are  given  in  Fig.  3.3  and  3.4  for  the  mass  mean  diameters.  The  mass  mean  diameters 
monitored  along  the  wall  and  axis  are  given  in  Fig.  3.3  and  values  measured  at  nozzle  exit  in  radial  direction  are  given 
in  Fig.  3.4. 

It  is  easily  seen  that  the  mass  mean  diameters  in  Fig.  3.3  are  slightly  sensitive  to  the  meshes  studied.  The  mass 
mean  diameter  varies  slightly  after  the  throat  following  the  axis  and  it  has  slightly  larger  value  on  265x80  than  132x60 
following  the  wall.  However,  the  mass  mean  diameter  differences  between  two  cases  are  approximately  less  than  three 
percent.  In  addition,  there  is  slightly  discernable  difference  of  mass  mean  diameters  at  nozzle  exit  between  two  cases. 
This  difference  is  observable  in  Fig.  3.4  at  the  center  and  wall  but  the  amount  of  difference  is  still  under  three  percent. 
It  looks  like  that  the  grid  should  be  chosen  carefully  at  around  the  center  line  after  the  throat  and  the  wall  for  more 
accurate  prediction.  However,  these  differences  do  not  result  in  the  averaged  value  of  particle  size  at  nozzle  exit  as  it 
is  shown  in  table  3.8,  because  the  differences  present  only  in  very  tiny  regions.  The  averaged  mass  mean  diameters  at 
nozzle  exit  plane  of  two  cases  are  0.878  and  0.879  pm. 


Table  3.8  Averaged  particle  characteristics  at  nozzle  exit  planes 


Grid  size . 

132x60 

265x80 

Dio 

0.213 

0.213 

D43 

0.878 

0.879 

Dm 

0.168 

0.168 

0.687 

0.687 

Figure  3.3  Axial  variation  of  the  mass  mean  diameter  along  the  wall  and  axis 


Radial  distance  from  the  axis,  y  (m) 

Figure  3.4  Radial  variation  of  the  mass  mean  diameter  at  nozzle  exit  planes 


3.3  Results  and  discussions 

3.3.1.  Comparison  between  the  prediction  and  Crowe  et  al.  [50]  ’s  experiments 

A  series  of  simulations  was  carried  out  to  predict  mass  mean  diameter  at  nozzle  exit  in  the  nozzle  configuration 
given  by  Crowe  et  al.[50].  As  described  above,  the  small  solid  rocket  motors  are  fired  into  the  collection  tank  with  the 
nozzle  [50]  and  without  nozzle  [51]  to  assess  the  particle  conditions  obtained  directly  from  the  chamber.  Ther  results 
are  presented  in  Figure  3.5.  The  chamber  pressure  was  changed  by  varying  the  propellant  burning  area  with  a  fixed 
nozzled  attached.  The  experimental  chamber  pressures  are  approximately  320,  570,  900,  and  1000  psi.  According  to 
Crowe  et  al.  [50],  they  observed  the  chamber  pressure  variation  during  test,  therefore,  we  chose  the  pressure  values 
from  Hermsen  [10]  to  account  for  the  pressure  variation.  The  chamber  pressures  of  simulations  are  470,  650,  and  980 
psi  and  the  corresponding  case  numbers  are  from  C-l  to  3,  respectively.  The  averaged  mass  mean  diameter  at  nozzle 
exit  plane  is  obtained  from  the  simulation.  The  predicted  results  are  presented  in  figure  3.5  and  we  observed  a  good 
agreement  between  the  measured  and  predicted  particle  size  over  500  psi  chamber  pressure.  In  case  of  650  psi 
chamber  pressure,  the  predicted  mass  mean  diameter  is  slightly  smaller  than  the  minum  experimental  value  at  570  psi 
chamber  pressure.  However,  the  mass  mean  diameter  used  in  the  simulation,  0.74  pm,  corresponds  to  the  lowest 
measured  mass  mean  diameter  directly  obtained  from  the  motor.  We  expect  to  obtain  the  larger  mass  mean  diameter 
using  the  maximum  measured  mass  mean  diameter  in  the  spread.  At  the  lower  pressure  level,  470  psi,  the  predicted 
mass  mean  diameter  is  approximately  30%  less  than  the  measured  value  at  320  psi.  Although  the  other  collision 
mechanisms  which  are  not  included  in  the  modeling  might  be  significant  in  this  range,  the  predicted  value  is  still  in  35% 
error  range  of  Hermsen’ s  correlation. 


Chamber  Pressure,  Pc  (psi) 

Figure  3.5  Predicted  and  measured  mass  mean  diameter  at  nozzle  exit  planes 
of  particles  obtained  from  Crowe  et  al.[50] 


3.3.2  The  effect  of  chamber  pressure 

The  current  modeling  described  previously  was  compared  here  to  Hermsen’s  correlation  and  the  dependence  of 
mass  mean  diameter  on  the  chamber  pressure  level  was  analyzed.  A  series  of  simulations  was  performed  to  compare 
the  results  to  the  correlation  under  the  nozzle  configuration  used  by  Shegal[53]  for  150psi  chamber  pressure.  The 
corresponding  case  numbers  are  from  P-1  to  6. 

Figures  3.6-9  show  the  effect  of  the  chamber  pressure  on  the  mass  mean  diameter  along  the  axis  and  wall,  and  at 
nozzle  exit  plane.  Figure  3.6  shows  that  the  largest  growth  rate  ocurrs  in  the  convergent  section,  specially,  the  most 
growth  occurs  in  short  region  in  front  of  the  nozzle  throat.  The  mass  mean  diameter  after  passing  the  throat  shows 
little  growth  in  divergent  section.  The  breakup  mechanism  due  to  velocity  lags  might  be  balanced  with  collision,  then 
the  mass  mean  diameter  approaches  a  certain  value  in  this  region. 

Figure  3.7  shows  the  variation  of  mass  mean  diameter  in  the  boundary  layer  and  it  shows  more  complicated 
variation  than  centerline  case.  It  is  observed  that  there  exist  three  peaks  for  all  cases.  The  first  peak  corresponds  to  the 
region  which  the  geometrical  compression  starts  from.  At  this  corner,  the  significant  shearing  motion  of  mean  flow 
occurs  due  to  the  recirculation  then  it  is  responsible  for  the  growth  of  particles.  This  finding  is  significant  since  it 
explains  the  importance  of  smooth  geometry.  The  second  peak  corresponds  to  the  maximum  mass  mean  diameter 
regions  formed  by  turbulence  motion.  Here  we  did  not  include  the  results  obtained  using  only  turbulent 
collision/breakup  mechanism,  the  maximum  mass  mean  diameter  occurs  at  this  region  as  we  reported  it  in  the 
previous  report.  The  interesting  results  about  this  peak,  the  collision  effects  due  to  turbulence  motion  are  less  effective 
with  larger  chamber  pressure  (or  collision  due  to  mean  flow  is  more  effective  than  turbulence  motion).  It  is  clearly 
observed  that  the  area  taken  by  this  peak  is  wider  with  low  pressure  than  high  pressure  case.  The  third  peak  is  the 
location  which  the  breakup  process  becomes  dominant  due  to  the  velocity  lag.  As  it  is  shown  in  the  figure,  the 
decrease  of  mass  mean  diameter  at  this  region  is  more  significant  in  high  chamber  pressure  case  than  low  pressure 
case.  The  reduction  rate  is  large  with  high  pressure  at  this  region.  The  interesting  observation  is  that  this  peak  is 
located  after  the  throat  with  low  chamber  pressure  and  located  before  throat  with  high  chamber  pressure.  Finally,  it  is 
observed  that  the  variation  of  mass  mean  diameter  in  the  short  region  at  around  nozzle  exit. 

Figure  3.8  shows  the  variation  of  the  mass  mean  diameter  in  radial  direction  at  nozzle  exit  plane.  It  is  clearly 
observed  that  the  large  amount  of  growth  occurs  within  boundary  layer  and  it  occurs  more  in  high  pressure  case  than 
low  pressure  case.  However,  the  mass  mean  diameter  is  almost  constant  in  most  region  in  this  direction.  The  growth 
occurring  within  center  region  is  mostly  due  to  the  mean  flow  shearing  as  the  velocity  in  radial  direction  increases 
radially. 

The  averaged  value  of  particle  size  at  nozzle  exit  is  summarized  in  table  3.9  and  the  predicted  mass  mean  diameter 
at  nozzle  exit  is  compared  with  Shegal’s  experiments,  Crowe  and  Willoughby[52]’s  prediction,  and  Hermsen’s 
correlation  in  figure  3.9.  The  predicted  results  are  mostly  less  than  the  measured  results  by  Shegal.  However,  Dobbins 
and  Strand  [55]  lately  indicated  that  Shegal’s  experimental  results  did  not  agree  with  other  measurements.  Dobbins 


and  strand  found  that  the  particle  size  increases  by  a  factor  of  1.7  with  a  ten-fold  increase  while  Shegal’s  experimental 
results  gave  increases  by  a  factor  of  5  with  a  ten-fold  increase,  approximately.  Therefore,  it  may  not  be  meaningful  to 
compare  the  prediction  with  Shegal’s  results.  Crowe  and  Willoughby’s  calculation  considering  the  slip  velocity 
between  the  particle  and  surrounding  and  the  momentum  exchange  in  collision  also  had  the  less  values  than  Shegal’s 
results. 

The  nozzle  inlet  conditions  reported  by  Fein[54]  for  Shegal’s  150  psi  case  are  0.240  of  the  particle  concentration, 
4.15  of  the  particle  residence  time,  and  0.79  pm  of  the  mass  mean  diameter.  Using  these  initial  conditions,  the 
simulation  performed  for  various  pressures  and  Hermsen’s  correlation  is  calculated.  Over  all  chamber  pressures,  the 
variation  trend  is  much  similar  to  Hermsen’s  correlation,  but  the  predicted  results  are  larger  than  the  results  from 
Hermsen’s  correlation.  However,  the  predicted  values  are  within  the  error  bounds  of  Hermsen’s  correlation  (35%) 
over  500  psi  chamber  pressure  and  slightly  over  35%  under  500  psi.  In  addition,  as  explained  above,  Shegal  reported 
the  larger  mass  mean  diameter  than  other  studies.  Because  he  used  the  same  technique  to  obtain  the  particle  size 
directly  from  the  motor,  his  results  might  have  a  larger  particle  size  than  the  actual  size.  Therefore,  we  expect  that  the 
resultant  particle  size  using  the  actual  inlet  particle  size  may  completely  fall  in  the  error  bounds  of  Hermsen’s 
correlation. 


Table  3.9  Averaged  particle  characteristics  at  nozzle  exit  planes 


Case  no. 

P-1 

P-2 

P-3 

P-4 

P-5 

P-6 

Dio 

0.218 

0.224 

0.232 

0.237 

0.243 

0.250 

D43 

0.971 

1.125 

1.299 

1.494 

1.693 

1.908 

Dm 

0.170 

0.172 

0.174 

0.175 

0.176 

0.178 

°s 

0.706 

0.733 

0.758 

0.783 

0.805 

0.824 

Axial  distance  from  the  throat,  x  (m) 

Figure  3.6  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  axis 


Figure  3.7  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  wall 
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Figure  3.8  Radial  variation  of  the  predicted  mass  mean  diameter  at  nozzle  exit  planes 
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Figure  3.9  Predicted  and  measured  mass  mean  diameter  at  nozzle  exit  planes  in 
Shegal[53]’s  experimental  condition  -  the  effect  of  chamber  pressure 

3.3.3  The  effect  of  particle  concentration  in  chamber 

The  current  modeling  described  previously  was  compared  here  to  Hermsen’ s  correlation  and  the  dependence  of 
the  particle  concentration  in  the  chamber  was  analyzed.  A  series  of  simulations  was  performed  to  compare  the  results 
to  the  correlation  under  the  nozzle  configuration  used  by  Shegal[53]  for  150psi  chamber  pressure.  Increasing  the 
Aluminum  loading  in  propellants  results  in  the  large  number  of  particles  in  chamber,  consequently,  it  has  more 
possibility  of  collision  passing  through  nozzle.  The  corresponding  case  numbers  are  from  M-l  to  6. 

Figures  3.10-13  show  the  effect  of  the  particle  concentration  on  the  mass  mean  diameter  along  the  axis  and  wall, 
and  at  nozzle  exit  plane.  All  of  figures  show  the  same  trends  obtained  in  the  chamber  pressure  variation  cases.  An 
interesting  observation  from  figure  3.11  is  that  the  first  peak  corresponding  to  the  region  of  the  corner  shows  very 
little  variation  according  to  the  variation  of  particle  concentration.  The  third  peak  is  located  almost  exactly  at  the 
throat. 

The  averaged  value  of  particle  size  at  nozzle  exit  is  summarized  in  table  3.10  and  the  predicted  mass  mean 
diameter  at  nozzle  exit  is  compared  with  Hermsen’ s  correlation  in  figure  3.13.  The  variation  trend  is  much  similar  to 
Hermsen’ s  correlation  over  all  particle  concentrations  used  in  simulations,  but  the  predicted  results  are  larger  than  the 
results  from  Hermsen’ s  correlation.  Below  the  particle  concentration  0.2,  the  predicted  values  are  slightly  larger  than 
35  %  of  Hermsen’s  correlation  and  the  predicted  values  are  approximately  35%  larger  than  Hermsen’s  correlation  for 
larger  particle  concentrations  than  0.2.  As  discussed  in  previous  section,  Shegal’s  results  might  have  a  larger  particle 
size  than  the  actual  size.  Therefore,  we  expect  that  the  resultant  particle  size  using  the  actual  inlet  particle  size  may  be 
within  the  error  bounds  of  Hermsen’s  correlation. 


Table  3.10  Averaged  particle  characteristics  at  nozzle  exit  planes 


Figure  3.10  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  axis 


Figure  3.11  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  wall 


Radial  distance  from  the  axis,  y  (m) 

Figure  3.12  Radial  variation  of  the  predicted  mass  mean  diameter  at  nozzle  exit  planes 


4.5 

4 

B 

3  3.5 


3 

o 

S 

a 

5 


a 

<D 


2.5  E- 
2 


1.5 


i  E- 


0.5 


: 

■ 

I  ♦  Crowe  and  Willough 

+  Current  Simulation 

>y’s  calculatioi 

7  1 

Current  Simulation  Crowe 

Pt.=550 
r„=4. 1 5 

and  Wilk 

Pt=600 
rf=15 
Dr 0.5 

tughby 

-r 

: 

D, =0.904 

Hermsen’s  corr.  fo 
Crowe  and  Willougt 

r 

lby 

- 

— 

-  - 

r 

/ 

-O' 

: 

** 

M  I 

fc  4 

◄ 

► 

◄ 

r 

j  ◄ 

r 

Herm 
-  curre 

sen  s  corr.  for  \ 

nt  simulation  ******** 

t 

1 

r 

m  4 

L- 

1 

__  Inlet  par 
7  current 

tide  size 
simulatio 

of 

n 

111 

u 

— 

particle  si 
and  Willc 

— 

i _ i _ I 

1II1C 

Crowe 

ze  of  — 
>ughby 

1  ■  ■  .  . 

0.1  0.15  0.2  0.25  0.3  0.35 

Particle  concentration  in  chamber,  cc  (g-mole/lOOg) 

Figure  3.13  Predicted  and  measured  mass  mean  diameter  at  nozzle  exit  planes  in 

Shegal[53]’s  experimental  condition  -  the  effect  of  particle  concentration 


3.3.4  The  effect  of  nozzle  scale 

A  last  series  of  simulations  were  performed  to  assess  the  effect  of  nozzle  scale  on  the  mass  mean  diameter  using 
the  nozzle  configuration  used  by  Shegal[53]  for  150psi  chamber  pressure.  As  discussed  in  Crowe  and  Willoughby [52], 
an  increase  in  nozzle  scale  gives  longer  particle  residence  time  in  the  nozzle,  which  implies  more  growth.  The 
corresponding  case  numbers  are  from  S-l  to  6. 

Figures  3.14-17  show  the  effect  of  the  nozzle  scale  on  the  mass  mean  diameter  along  the  axis  and  wall,  and  at 
nozzle  exit  plane.  All  of  figures  show  the  similar  trends  obtained  in  the  chamber  pressure  variation  cases.  An 
interesting  observation  from  figure  3.15  is  that  the  second  peak  of  S-6  shows  very  large  decrease  of  mass  mean 
diameter  in  boundary  layer.  The  second  peak  becomes  more  noticeable  as  the  nozzle  scale  increases.  From  figure  3.16, 
it  is  observed  that  the  large  amount  of  growth  occurs  within  boundary  layer  and  it  occurs  more  in  large  scale  nozzle 
than  small  scale. 

The  averaged  value  of  particle  size  at  nozzle  exit  is  summarized  in  table  3.11  and  the  predicted  mass  mean 
diameter  at  nozzle  exit  is  compared  with  Hermsen’s  correlation  in  figure  3.17.  The  variation  trend  is  very  similar  to 
Hermsen’s  correlation  in  small  scale  nozzles,  as  the  mass  mean  diameter  increases  with  an  increase  of  nozzle  scale. 
However,  it  is  observed  that  the  mass  mean  diameter  decreases  slightly  with  an  increase  of  nozzle  scale  in  large  scale 
nozzles.  Over  all  nozzle  scales  except  the  mimum  of  scales,  the  results  are  in  35%  error  bounds  of  Hermsen’s 
correlation. 


Table  3.11  Averaged  particle  characteristics  at  nozzle  exit  planes 


Case  no. 

S-l 

S-2 

S-3 

S-4 

S-5 

S-6 

Dio 

0.243 

0.255 

0.273 

0.284 

0.286 

2.830 

D43 

1.693 

1.943 

2.261 

2.419 

2.474 

2.437 

Dm 

0.176 

0.182 

0.192 

0.198 

0.199 

0.198 

°s 

0.805 

0.823 

0.840 

0.845 

0.848 

0.847 

Normalized  axial  distance  from  the  throat,  x/Dt 
Figure  3.14  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  axis 
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Normalized  axial  distance  from  the  throat,  x/Dt 
Figure  3.12  Axial  variation  of  the  predicted  mass  mean  diameter  along  the  wall 


Nonnalized  radial  distance  from  the  axis ,y/Dt 

Figure  3.13  Radial  variation  of  the  predicted  mass  mean  diameter  at  nozzle  exit  planes 


Nozzle  throat  diameter,  Dt  (m) 

Figure  3.14  Predicted  and  measured  mass  mean  diameter  at  nozzle  exit  planes  in 
Shegal[53]’s  experimental  condition  -  the  effect  of  nozzle  scale 


Conclusions 


A  new  model  is  under  development  to  assess  coalescence  and  breakup  processes  in  rocket  combustor  and  nozzle 
environments.  The  one-way  coupled  population  balance  equation  describing  the  change  of  number  concentration  by 
the  modeled  particle  to  particle  interactions  and  aerodynamic  forces  is  solved  using  the  direct  quadrature  method  of 
moments  (DQMOM).  The  required  parameters  to  describe  the  collision  and  breakup  processes  are  modeled  in  laminar 
and  turbulent  flow. 

The  modeling  was  compared  to  experiments  and  correlation  with  respects  to  the  variations  in  chamber  pressure, 
particle  concentration  in  chamber,  and  nozzle  scale.  The  comparisons  show  that  the  predicted  mass  mean  diameters 
are  in  a  good  agreement  with  experiments  and  correlation  over  500  psi  chamber  pressure.  The  predicted  mass  mean 
diameters  also  have  a  good  agreement  with  correlation  over  0.2  mole/100  g  particle  concentration  and  within  all  tested 
nozzle  scales.  These  results  indicate  the  validity  of  the  current  model  for  particle  growth/reduction. 

Coalescence  is  shown  to  occur  in  the  convergent  section  leading  to  the  throat,  while  breakup  processes  tend  to 
become  important  in  the  throat  region  and  exit  cone.  In  addition,  the  modeling  shows  that  more  growth  occurs  in 
boundary  layers  than  mean  flow  regions. 

The  restriction  of  current  model  is  the  necessity  of  the  accurate  information  on  the  particle  characteristics  in  the 
chamber.  Therefore,  the  analytical  techniques  accurately  addressing  the  particle  size  and  concentration  in  chamber 
shoud  be  obtained  for  more  accurate  predictions. 
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On  the  Dynamic  Response  of  Rocket  Swirl  Injectors 
Part  I.  Theoretical  Description  of  Wave  Reflection 

and  Resonance 


Maksud  (Max)  Ismailov1  and  Stephen  D.  Heister2 
Purdue  University.  701  W.  Stadium  Ave.,  West  Lafayette,  IN  47907 


Linear  analyses  are  developed  to  investigate  theoretically  how  disturbance  waves  are 
reflected  and  transmitted  in  the  vortex  chamber  of  a  classical  swirl  injector.  The  de¬ 
pendence  of  the  magnitude  of  the  wave  reflection  process  on  the  disturbance  frequency 
is  derived,  and  it  is  shown  that  this  dependence  may  exhibit  distinct  maximum  values. 
It  is  explained  that  the  frequencies  at  which  maximum  response  occurs  are  termed  the 
resonant  frequencies  of  the  swirl  injector.  In  general,  resonant  conditions  will  depend 
not  only  on  the  geometry  of  the  injector,  but  also  on  the  particular  flow  conditions. 
In  other  words,  for  a  given  injector  geometry,  there  are  specific  flow  conditions  that 
may  produce  resonance.  A  simple  formula  is  derived  for  the  primary  resonance  which 
corresponds  to  a  1/4  wave  oscillation  within  the  vortex  chamber.  Two  different  reso¬ 
nance  theories  are  presented,  which  vary  in  their  level  of  accuracy  of  description  of  the 
flow  transition  from  the  vortex  chamber  to  the  nozzle.  Results  are  provided  for  both 
of  these  models. 
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c 


angular  momentum  constant  ( C  =  u0r) 


f  =  disturbance  frequency 


kn  =  wave  number  in  uniform  nozzle  region 

Lc  =  length  of  conical  convergence  section 

Ln  =  length  of  nozzle 

Lv  =  length  of  vortex  chamber 

Nin  =  number  of  tangential  inlets 

R  =  radius  of  solid  boundary 

Rin  =  inflow  radius  ( Rin  =  Rv  —  Rt) 

Rn  =  radius  of  nozzle 

Rt  =  radius  of  tangential  inlet 

Rv  =  radius  of  vortex  chamber 

The  =  steady  free  surface  radius  at  head  end,  z  —  0 

r„  =  steady  free  surface  radius  in  uniform  nozzle  region 

tv  =  steady  free  surface  radius  in  uniform  vortex  chamber  region 

u0  =  circumferential  velocity 

uz  =  axial  velocity 

uzn  =  axial  velocity  in  uniform  nozzle  region 

uzv  =  axial  velocity  in  uniform  vortex  chamber  region 

Win  =  tangential  inlet  inflow  velocity 

a  =  angle  of  solid  wall  convergence 

8  =  steady  free  surface  radius  at  any  axial  position  of  internal  flow 

tj  =  free  surface  deflection  away  from  8 

A  =  wave  length 

IT inj  =  total  injector  response 

u  =  angular  disturbance  frequency 

prime  ()r  =  fluctuation  value  of  parameter 

bar  ()  =  steady  state  value  of  parameter 

hat  ()  =  amplitude  of  fluctuation  of  disturbed  parameter 

star  ()*  =  dimensional  value  of  parameter 
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I.  Introduction 


Swirl  injectors,  or  simplex  atomizers,  are  used  in  a  variety  of  applications  ranging  from  agricul¬ 
tural  sprays  to  use  within  high  flow  aerospace  combustors.  In  particular,  these  injectors  have  seen 
wide  application  in  Russian  rocket  engines  and  in  numerous  gas  turbine  engines.  In  these  applica¬ 
tions,  it  is  well  known  that  the  injector  can  play  a  role  in  stabilizing  or  destabilizing  combustion 
processes  as  it  acts  as  an  active  element  within  the  system;  in  particular  the  injector  response  at 
acoustic  modes  of  the  chamber  is  highly  critical.  Should  the  injector  exhibit  a  resonance  condition 
at  a  fundamental  chamber  acoustic  frequency,  the  potential  exists  for  rapid  growth  of  combus¬ 
tion  instabilities.  For  this  reason,  the  dynamic  response  of  swirl  injectors  is  a  topic  of  significant 
importance  in  this  community. 

Despite  this  interest,  there  are  rather  limited  works  that  discuss  the  dynamic  response  of  swirl 
injectors.  The  most  widely  recognized  analytical  work  in  this  field  is  the  1979  book  by  Bazarov  [1], 
who  developed  a  linear  theory  based  on  small  disturbances  propagating  through  the  injector  flow. 
Overall,  Bazarov’s  theory  [1,  2]  presents  itself  as  a  valuable  analytical  tool  that  may  indicate  on 
which  frequencies  what  strength  of  injector  response  one  can  expect.  Before  Bazarov’s  work,  there 
did  not  exist  any  theory  similar  to  this  one  in  its  level  of  extensive  description  of  injector  dynamics. 
Because  this  is  the  current  standard  for  use  in  analysis  of  dynamic  response  of  swirl  injectors,  we 
developed  an  in-house  code  to  be  able  to  use  this  as  a  predictive  capability  within  our  group.  The 
code  was  validated  against  results  published  within  Bazarov’s  book.  The  code  was  exercised  over 
a  range  of  conditions  in  order  to  assess  its  behavior  [3,  Chap.  2].  As  a  result  of  these  efforts,  we 
identified  several  aspects  of  the  theory  which  provide  for  motivation  for  improvements  and  further 
study: 

Regarding  surface  wave  treatment: 

•  The  variation  of  internal  flow  boundaries  is  simplified  to  the  case  of  sudden  film  thickness 
change,  as  the  vortex  chamber  transitions  to  the  nozzle,  which  permits  one  to  ignore  wave 
refraction  and  results  in  a  reflection  coefficient  that  does  not  depend  on  the  disturbance  fre¬ 
quency. 

•  When  treating  surface  waves,  that  are  reflecting  back  an  forth  in  the  vortex  chamber,  the  wave 
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speed  taken  for  both  downstream  and  upstream  traveling  disturbances  is  such  that  it  is  valid 
for  downstream  traveling  disturbances  only. 

•  Only  long  wave  speed  relations  are  considered,  that  do  not  formally  apply  at  high  disturbance 
frequencies,  where  Kelvin’s  [4]  dispersion  becomes  more  appropriate. 

Regarding  vorticity  wave  treatment: 

•  The  process  of  liquid  issuing  from  the  injector  tangential  inlet  into  the  vortex  chamber  is 
represented  by  a  conformal  mapping  of  the  half-strip  into  the  half-plane,  however,  no  clear 
theoretical  evidence  exists  justifying  this  representation. 

•  To  determine  the  phase  shift  in  the  radial  direction,  the  expression  for  steady  state  radial 
velocity  is  used  which  does  not  obey  Laplace’s  equation  describing  the  potential  flow  at  the 
steady  state. 

•  To  compute  the  time  lag  in  the  radial  direction,  the  radial  distance  from  the  cylindrical  wall 
of  the  vortex  chamber  to  the  point  of  interest  at  some  arbitrary  radius  is  divided  by  the 
steady  radial  velocity  corresponding  to  this  point,  however,  it  is  not  taken  into  account  that 
this  radial  velocity  actually  varies  with  radius,  so  that  an  integral  expression  for  the  time  lag 
should  be  used. 

These  issues  provided  a  motivation  for  the  present  study. 

More  recently,  Park  [5]  investigated  the  dynamics  of  swirl  injectors  by  using  axisymmetric 
boundary  element  method  (BEM),  which  assumes  that  the  How  is  incompressible  and  irrotational. 
The  way  Park  modeled  the  How  unsteadiness  was  by  fluctuating  sinusoidally  the  inflow  velocity, 
while  the  pressure  inside  the  core  remained  constant.  His  work  was  among  the  first  that  attempted 
to  model  the  unsteady  dynamics  of  swirl  injector  by  methods  of  CFD.  Richardson  [6,  7]  continued 
the  above  work  of  Park  [5]  and  relaxed  the  condition  of  constant  gas  pressure  in  the  core  and  allowed 
it  to  fluctuate.  Now,  the  inflow  velocity  was  computed  based  on  instantaneous  pressure  drop  and 
the  radius  of  the  core  surface  at  the  head  end.  Recently,  a  research  group  in  Korea,  Khil  et  al.  [8],  [9] 
pulsed  the  flow  in  the  range  of  frequencies  of  100-300  Hz,  and  presented  the  experimental  data  for 
the  pressures  and  mass  flow  rates.  Experimental  results  have  been  also  published  by  Ahn  [10,  11], 
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whose  work  is  aimed  at  achieving  response  data  in  the  broader  frequency  range,  possibly  the  same 
as  in  Bazarov’s  [1]  experiments. 

Figure  1  provides  a  schematic  of  a  classical  swirl  injector.  An  “air  core’*  evolves  at  the  center 
of  the  injector  due  to  the  vortical  inflow.  The  shape  of  the  free  surface  shown  schematically  in 
Fig.  1  has  a  head-end  transition  region  that  results  from  the  flow  accelerating  in  the  axial  direction 
from  its  initial  tangential  injection.  A  similar  transition  region  exists  near  the  nozzle  as  the  flow  is 
accelerated  via  the  conical  convergence  section.  One  of  the  major  hurdles  one  must  undertake  to 
assess  the  dynamic  response  of  a  swirl  injector  relates  to  the  propagation,  refraction,  and  reflection 
of  waves  that  form  on  the  free  surface  of  the  liquid  film.  Complex  dynamic  patterns  can  result  from 
interaction  of  these  waves,  and,  as  the  flow  is  swirling,  the  undulations  can  also  lead  to  local  changes 
in  the  swirl  level.  For  these  reasons,  the  analysis  of  the  dynamic  response  of  these  injectors  is  quite 
challenging.  In  this  study,  we  focus  on  the  disturbance  wave  reflections  and  the  resonance  caused 
by  them.  There  is  no  indication  in  the  existing  literature  that  this  topic  has  been  investigated  to 
date.  Note  that,  at  the  present  time,  all  analytical  models  describing  the  steady-state  injector  free 
surface  do  not  consider  the  head  end  and  the  nozzle  entrance  transition  regions  at  all,  because  the 
engineering  design  interest  lies  mainly  in  finding  the  core  radii  in  the  uniform  regions  of  the  injector 
(Taylor  [12],  Bayvel  and  Orzechowski  [13],  Bazarov  [1]).  In  reality,  the  variation  of  the  free  surface 
shape  in  the  said  transition  regions  does  take  place,  and  is  smooth  and  continuous,  which  sets  the 
favorable  conditions  for  the  wave  refractions  and  reflections  to  occur  [3]  as  follows  next. 

Imagine  that  an  incident  wave  has  been  induced  in  the  vortex  chamber  and  is  moving  in  the 
positive  axial  direction  towards  the  nozzle.  At  the  entry  to  the  conical  convergence  section,  an 
abrupt  change  in  the  bulk  flow  velocity  and  the  liquid  film  thickness  takes  place.  Thus,  when  an 
incident  wave  arrives  at  that  discontinuity  plane,  one  portion  of  it  reflects  back  into  the  vortex 
chamber,  in  the  negative  axial  direction,  and  another  portion  transmits  further  into  the  transition 
region.  Now,  let  us  track  the  reflected  portion  of  it  in  the  vortex  chamber.  When  the  reflected  wave 
comes  to  the  head  end  solid  wall,  it  is  reflected  once  again,  in  the  positive  axial  direction.  Thus, 
we  see  two  planes  generating  a  wave  reflection,  suggesting  that  a  standing  wave  pattern  may  be 
arranged  in  the  vortex  chamber. 
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Fig.  1  Swirl  injector  flow  schematic 


The  amplitude  of  the  resulting  standing  wave  depends  on  the  reflection  and  transmission  char¬ 
acteristics  quantifying  the  events  of  reflections  just  described.  Because  these  characteristics  depend 
on  the  disturbance  frequency,  then,  there  may  be  frequencies,  at  which  this  amplitude  is  maximized. 
These  frequencies  are  the  natural  frequencies  of  the  swirl  injector.  Since  the  swirl  injector  is  coupled 
to  the  combustion  chamber  of  the  rocket  engine,  then  a  resonance  may  occur  when  the  combustion 
instability  frequency  coincides  with  one  of  the  natural  frequencies  of  the  injector.  Accordingly,  fur¬ 
ther  on,  we  will  refer  to  the  natural  frequencies  of  the  swirl  injector  as  the  resonant  frequencies. 
For  this  reason,  it  is  highly  desirable  to  be  able  to  predict  these  resonant  conditions  and  that  is  the 
prime  motivation  for  the  present  study. 

We  can  treat  this  problem  by  using  three  different  levels  of  approximation,  depending  on  what 
parts  of  the  injector  are  included  in  the  analysis.  In  the  first  approximation,  Sec.  IV,  we  consider  the 
vortex  chamber  and  the  nozzle  by  assuming  that  they  are  connected  by  a  sharp  (radial)  step,  and 
look  for  the  nozzle  effect  on  the  wave  reflection  and  resonance  characteristics.  We  shall  refer  to  this 
treatment  as  the  Abrupt  Convergence  Resonance  Model,  or  ACRM.  In  the  second  approximation, 
Sec.  V,  we  consider  all  three  elements  of  the  swirl  injector,  the  vortex  chamber,  the  conical  convergent 
section,  and  the  nozzle,  and  look  for  their  collective  influence  on  the  wave  reflection  and  resonance. 
We  shall  refer  to  this  as  the  Conical  Convergence  Resonance  Model,  or  CCRM. 

In  both  approximations,  we  approach  this  problem  by  matching  the  instantaneous  mass  flow 
rate  and  momentum  fluctuations  at  locations  where  the  flow  has  discontinuities.  For  this  purpose, 
in  Sec.  Ill,  the  expressions  for  the  instantaneous  mass  flow  rate  and  momentum  fluctuation  in  terms 
of  the  disturbance  frequency  will  be  derived. 

Both  of  these  approximations  are  equally  important,  because  each  one  of  them  may  be  used  for 
the  assessment  of  the  swirl  injector  in  terms  of  its  resonant  characteristics,  depending  on  the  level 
of  detail  required  in  the  assessment.  Also,  they  may  serve  as  a  cross  check  between  each  other,  since 
this  topic  is  new  and  the  results  have  to  be  anchored  to  one  another  when  possible. 

Note  that  in  all  models  considered,  we  will  assume  that  the  circumferential  velocity  follows  the 
free  vortex  law  as  ug  =  C/r.  Finally,  we  shall  emphasize  that  we  will  deal  with  the  long  waves  only. 


II.  Fundamental  Condition  for  Resonance  from  Wave  Considerations 


We  know  that  a  general  swirl  injector  with  an  arbitrary  angle  of  conical  convergence  generates 
an  oblique  reflection  at  the  point  where  the  vortex  chamber  connects  to  the  conical  section.  For 


simplicity,  let  us  imagine  that  we  have  placed  a  straight  step  discontinuity  at  that  point,  so  that 
all  reflections  in  the  vortex  chamber  become  normal.  Further,  let  us  completely  disregard  the  part 
of  the  injector  downstream  of  that  step  discontinuity.  Figure  2  then  shows  the  assumed  injector 
representation  for  this  problem. 
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Fig.  2  Schematic  of  injector  representation  for  ACRM 


Now,  let  us  imagine  what  would  the  standing  wave  pattern  look  like  when  the  injector  resonates 
with  downstream  processes  (such  as  combustion).  Consider  vortex  chamber  resonance  from  the 
injector  response  perspective.  First,  we  expect  that,  at  resonance,  the  injector  response  should  be 


at  its  maximum.  Second,  from  its  definition  (Bazarov  [1]), 


mn 


1ImJ  "  AP|n, 


(1) 


APi  nj 

where  ?7in  is  the  nozzle  mass  flow  rate  and  A pinj  is  the  injector  pressure  drop,  we  see  that  the 
injector  response  is  maximized  when  the  magnitude  of  the  nozzle  mass  flow  rate  fluctuation,  mj,, 
is  at  its  maximum,  and  the  magnitude  of  the  injector  pressure  drop  fluctuation,  A p[nj,  is  at  its 
minimum.  This  is  true,  if  and  only  if  we  have  a  node  at  the  head  end  and  an  antinode  at  the  nozzle 
entrance  [3,  Sec.  5.2].  Figure  3  shows  schematically  different  possible  modes  of  the  wave  pattern 
that  we  may  anticipate  when  the  injector  resonates. 
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Fig.  3  Schematic  of  standing  wave  pattern  when  swirl  injector  is  at  resonance 

In  addition  to  the  assumed  radial  step  discontinuity,  we  also  assume: 

1.  The  disturbance  wave  speed  is  much  larger  than  the  axial  bulk  flow  velocity  in  the  vortex 

chamber,  so  that  we  can  think  of  the  fluid  in  the  vortex  chamber  as  being  quiescent  with  zero 
_  duzv 

axial  velocity,  or  =  0  and  — - —  =  0. 

dz 

2.  Since  the  variation  of  the  free  surface  radius  at  the  head  end  is  small  [3,  Chap.  3],  we  can 
overall  assume  that  The  =  Tv  (see  Fig.  1). 
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In  accord  with  these  assumptions,  we  can  write  the  wave  equation  [1,  14],  which  describes  the 


height  77  of  long  waves  traveling  on  the  core  of  swirling  fluid  in  the  vortex  chamber,  as 

d2jl  _C2  RZ~Tv  d2f) 

dt2  It*  dz 2 


(2) 


wThere  the  injector  dimensions  are  shown  in  Figs.  1  and  3,  and  C  =  ugr  is  the  angular  momentum 
constant  in  the  swirling  flow.  The  general  bounded  solution  to  this  equation  is  a  periodic  function 

77  =  [P  cos  (k^t )  H-  Q  sin  (k^t )]  [A  cos  ( kz )  +  B  sin  ( kz )] 

o  _  r2 

where  P,  Q,  A,  and  B  are  the  unknown  constants,  j2  =  C~ — vn  v ,  and  k  is  the  wave  number. 

2r4 

The  positions  of  nodes  and  antinodes  in  Fig.  3  indicate  two  boundary  conditions  for  77  valid  at 
all  times: 

(a)  1 77I  should  be  zero  at  the  head  end,  z  =  0. 

(b)  1 77I  should  be  maximum  at  the  end  of  the  vortex  chamber,  z  =  Lv. 

From  (a),  we  conclude  that  >1  =  0,  which  reduces  the  general  solution  to 

77  =  [P  cos  (kjt)  +  Q  sin  (kjt )]  sin  (kz) 

where  we  have  absorbed  B  into  P  and  Q.  Now  the  standing  wave  pattern  is  clearly  seen.  From  (b), 
we  deduce  that  sin  ( kLv )  =  ±1,  which  yields  the  resonant  wave  numbers: 


k  =  77——,  n  =  l,  3,  5, . 
21/  ’  ’ 


(3) 


The  resonant  frequencies,  u.’0  =  k'y,  are  respectively 


n= w,.. 


where  by  subscript  zero  we  emphasize  the  notion  of  resonance. 


(4) 


III.  Long  Wave  Fluctuations  of  Mass  Flow  and  Momentum  in  Cylindrical  Flow  Sections 

In  this  section,  based  on  the  definitions  of  instantaneous  mass  flow  rate  and  momentum,  we  will 
derive  expressions  for  the  fluctuating  parts  of  each  of  them,  written  up  to  the  first  order  in  77,  that  are 
valid  in  a  purely  cylindrical  section  of  the  swurling  How,  in  which  the  radii  of  steady  How  boundaries 
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and  bulk  stream  velocity  are  considered  constant,  and  which  is  experiencing  long  wave  perturbations. 
We  need  to  know  this  information  for  use  in  the  more  precise  wave  reflection/resonance  models 
(ACRM  and  CCRM)  of  the  swirl  injector  that  will  come  later  in  the  next  two  sections. 

Let  us  assume  that  we  can  split  the  axial  flow  velocity,  uz,  and  the  cross-sectional  area  of  the 
gaseous  core,  A,  into  the  mean  and  disturbed  parts  as  uz  =  uz  +  u'z  and  A  =  A  +  A'  respectively, 
where  the  disturbances  are  assumed  small,  and  given  by  Fourier  waves 

u'z  =  uze^kz~ut)  A'  =  AeiV‘z~ut'>  T]  =  f,ei(kz~ut) 

Let  us  determine  uz  and  A1  in  terms  of  77.  From  the  definition  of  A, 

A  =  A  +  A'  =  7 r(<5  +  7])2  ~  nS2  +  27v5rj 


follows  that  A f  =  2ttSj].  To  find  uzf  the  cylindrical  constant  velocity  How  setup  assumed  here  in 
permits  us  to  use  the  linearized  continuity  equation  by  Darmofal  et  al.  [15],  which  in  combination 
with  the  expression  for  A'  gives: 


/  o  c  ku*  ~  " 

“•  =  um=e?i 

where  S  is  the  steady  free  surface  radius,  and  R  is  the  solid  boundary  radius. 


(5) 


A.  Long  Wave  Fluctuation  of  Mass  Flow  Rate 

By  definition,  the  instantaneous  mass  flow  rate  at  any  flow  cross-section  can  be  written  as 

R  r  R  5+r,  R  S+t}  “I 

I  /  Ti  Q/r rr/Jr  —  ^  c)'trrRr  -l_  /  9<r rrrlr  —  /  9 'rrr rl t  ! 

& 

Recognizing  the  steady  mass  flow  rate  in  the  first  integral,  neglecting  the  last  integral  as  it  is  a 
higher  order  term,  and  eliminating  other  higher  order  terms,  the  unsteady  part  of  the  mass  flow 
rate  becomes: 


H  R  R  d+T7 

h  =  p  j  (uz  +  ^)27rrdr  =  p  J  uz27rrdr  —  J  uz27rrdr  +  J  uz2nTdr  —  J  u'z2nrdr 
<5+17  La  5  s  s 


1 

j 


TO 7  =  —2'Kp8uzr\  +  pTT  ( R 2  —  52)  vlz 

We  can  modify  further  this  equation  by  substituting  the  velocity  fluctuation  from  Eq.  (5), 

ml  =  —2npS^-r]  (6) 
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B.  Long  Wave  Fluctuation  of  Momentum 

Assuming  there  are  no  external  body  or  friction  forces  acting  on  the  flow,  we  can  define  the  total 
momentum  of  the  flow,  Me  =  Me  +  A/£,  as  being  consistent  of  the  kinetic  part,  A/u  =  Mu  +  A/', 

and  the  pressure  part,  Mp  =  Mp  +  A/',  which  we  may  write  as 

R  R 

A/e  =  Mu  +  Mp  =  /?  /  (w2  +  )227rrdr  +  p(t)  2nrdr  (7) 

Let  us  investigate  both  of  these  integrals  separately.  Starting  from  the  kinetic  part,  we  may  write 
R 

Mu  —  A/u  +  A/'  =  p  I  (uz  +  u,z)22nrdr  =  pu2TT  ( R 2  —  S 2)  —  pu2n28r]  +  p2uzu'z7r  ( /? 2  —  <52) 

*+r? 

From  here  we  can  extract  the  unsteady  part,  and  then  substitute  the  velocity  fluctuation  from 
Eq.  (5)  as  follows: 


Mfu  =  —p2Tr5u2zj]  +  p2n  ( R 2  —  82)  uzv!z  =  pn  (  —2 8u2  +  4 uz8 


kuz  —  LJ 


Next,  since  the  flow  is  purely  axial  and  we  do  not  have  the  radial  velocity  fluctuations  in  the  long 
wave  treatment,  the  pressure  at  some  radius  in  the  How  may  be  written  as 


where  t  is  a  dummy  radius.  Now,  we  can  insert  this  expression  into  the  pressure  part  of  momentum 
Eq.  (7)  to  get 


A/p  =  A/p  +  A/'  = 


rl 

p  =  j  P{r) 


Inrdr  =  pnC2  j  [*2  "  (*+  V)2]  ~  1 


To  work  on  the  first  term  in  this  expression,  we  will  use  the  binomial  expansion 


1  1  2 


5 2  8* 


(<*+»?)■ 


By  substituting  into  previous  equation,  we  can  write 


2 


2  [R2  ~  (S  +  vf 


,1  R2-S2  R2-52  1 


25-  53  n  51] 


To  modify  the  logarithm  term  in  Eq.  (9),  we  will  use  the  binomial  expansion 


1  1  _  1 
J+^~5  S2*1 
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and  the  series  expansion  [16,  p.  Ill] 


/  \  (x  ~  l)2  (x  ~  l)3 

lnx  -  (x  -  1) - jt—  +  “ — 5 - 


which  we  write  up  to  the  third  term  in  this  study.  Note  that  this  expansion  is  formally  valid  for 
|x  -  1|  <  1,  or  0  <  x  <  2.  In  this  study,  for  simplicity,  we  wTill  assume  that  this  series  expansion  is 
valid  for  the  whole  range  of  the  argument,  or  for  x  >  0.  Ideally,  one  would  write  as  many  terms  in 
this  expansion  as  possible  to  make  the  end  result  more  precise.  Then,  the  logarithmic  term  can  be 
written  as: 


,  R  (R  \  1  (R  \2  1  (R  \3  R  (  Rr  0i?  \ 

n<S+r/  VJ  V  2(5  V  +3\J  V  +  <PV  52+3(S  7 


(11) 


Plugging  Eqs.  (10)  and  (11)  back  into  the  pressure  part  of  momentum,  Eq.  (9),  we  have 
Mp  —  pnC 2 


oq 

CN 

*-0 

1 

cs 

oq 

1 - 

->K( 

R-i 

y_ 

2S2  V  S 

5  1 

) 

sU  vj 

...  +  pnC2 

'  R2-S 2 

1 

-( 

'  R2  „R 

- - 1-  3  —  - 

52  5 

[  5 3 

"  S  - 

<52  \ 

>)] 


From  here,  we  conclude  that  the  unsteady  part  of  the  pressure  momentum  is  given  by 
R2-5 2 


M;  =  pnC2 


S3 


-  \  -  i  (~w + 3f  - 3)] " = pnC2¥ {R3  - AR2S + ™52)  11  (12) 


Adding  the  unsteady  kinetic  and  pressure  parts  given  by  Eqs.  (8)  and  (12),  and  rearranging,  we 
arrive  at  the  expression  for  the  total  momentum  fluctuation: 


M £  =  pTT 


-2 Sul  +  4g,^fctt*.  -  +  C2±  (R3  -  4 R2S  +  3 R52) 

K  O'* 


(13) 


IV.  Abrupt  Convergence  Resonance  Model  (ACRM):  Wave  Reflections  and  Resonance 
when  Vortex  Chamber  and  Nozzle  are  Connected  with  an  Abrupt  Step  Discontinuity 

In  this  section,  let  us  build  upon  the  above  mentioned  idea  that  we  have  an  abrupt  step  change 
at  the  end  of  the  vortex  chamber,  and  connect  the  nozzle  just  downstream  of  that  step,  with  its 
own  solid  boundary  and  free  surface  radii,  Rn  and  rn  respectively,  Fig.  4.  Consider  an  incident 
downstream  traveling  wave  D  exp  [i  (k2z  —  ut)]  in  the  vortex  chamber.  When  it  comes  to  the  step 
discontinuity,  a  part  of  it  reflects  back  as  an  upstream  traveling  wave  B  exp  [i  (kjz  —  u-'f)]>  an^ 
another  part  of  it  propagates  further  into  the  nozzle  as  a  transmitted  downstream  traveling  wave 
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Fig.  4  Schematic  of  wave  reflection  and  transmission  for  ACRM-2 


Aexp[i(fcn2  —  c tit)].  The  arrows  in  Fig.  4  indicate  the  respective  directions  in  which  all  of  these 
waves  travel,  and  k  1?  k2 ,  and  kn  denote  the  corresponding  wave  numbers.  Here  and  further  on,  the 
upstream  traveling  waves  will  be  indicated  by  odd  numbers,  and  the  downstream  traveling  waves 
-  by  even  numbers.  Note  that  there  may  not  exist  an  upstream  traveling  wave  in  the  nozzle  as 
the  local  How  speed  typically  exceeds  the  upstream  wave  speed.  Finally,  similarly  to  Section  II,  we 
will  assume  again  that  the  bulk  flow  axial  velocity  in  the  vortex  chamber  is  negligibly  smaller  in 
comparison  with  the  disturbance  wave  speed,  and  may  be  set  to  zero,  uzv  =  0. 

Based  on  Eq.  (5),  at  any  point  in  the  vortex  chamber,  we  may  write  the  collective  velocity 
fluctuation  of  the  waves  B  and  £>,  that  are  superposed  on  each  other,  as 


k\uzv  -  Lu 


k2uzv  —  uj 

*2  (J?l  -  r5) 


Qei(k2z-ut) 


At  the  head  end  solid  wall,  z  =  0,  we  know  that  the  axial  velocity  fluctuation  should  be  zero  at  all 
times.  Then,  noting  that  uzv  =  0,  we  can  reduce  this  equation  to 


u'z  (°. *)  =  -j^  +  =  0 


Again,  since  uzv  =  0,  the  wave  numbers  k \  and  k2  are  symmetrical,  k\  =  —k2.  But,  based  on  the 
last  equation,  this  means  that  D  =  B.  Conclusively,  we  can  say  that  we  have  a  purely  standing 
wave  in  the  vortex  chamber,  which  excites  an  outgoing  wave  of  amplitude  A  in  the  nozzle. 
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To  relate  A,  to  the  amplitude  of  the  standing  wave  in  the  vortex  chamber,  D ,  consider  the 
balance  of  the  mass  How  rate  fluctuation  at  the  matching  boundary,  z  =  Lv  (see  Fig.  4).  Based  on 
the  derived  Eq.  (6)  for  the  mass  How  rate  fluctuation,  we  may  write 

—2nprv  \^-Be^k'L~-^  +  =  -2n prn^-Aei(k"L“~wt) 

k i  ko  kn 

By  substituting  the  equalities  D  =  B  and  k i  =  — k 2,  and  rearranging,  we  can  obtain  the  expression 
for  A  in  terms  of  D : 

A  =  ±2 - - D  (14) 

-LJLPiknLv 

K 

What  does  this  equation  mean  from  the  injector  resonance  point  of  view?  In  this  equation, 
k2  and  kn  depend  on  the  disturbance  frequency  through  the  long  wave  speed  relationship,  given  in 
general  as 

u-ku,  =  ±k][c^^-  (15) 
and  which  for  the  vortex  chamber  and  the  nozzle  is  given  by 

u=+k2)r^fi  u-knuzn  =  +kn 

These  equations  take  into  account  the  respective  film  thicknesses  in  the  vortex  chamber,  Rv  —  rv , 
and  the  nozzle,  Rn  —  rn.  Hence,  if  D  is  fixed  in  Eq.  (14),  the  magnitude  of  .4  should  vary  with 
regard  to  cj.  Then,  there  may  be  frequencies,  where  \A\  is  maximized,  thereby  causing  the  most 
pronounced  mass  flow  rate  pulsation  in  the  nozzle.  Theoretically,  these  frequencies  should  coincide 
with  the  resonant  frequencies,  Eq.  (4),  which  we  have  derived  in  the  fundamental  condition  for 
resonance  above  (Sec.  II). 

V.  Conical  Convergence  Resonance  Model  (CCRM):  Wave  Reflections  and  Resonance  in 
Injector  with  a  Conical  Convergence  Section 

Let  us  add  some  of  the  more  realistic  features  to  the  previous  model  by  considering  that  there 
are  the  following  additional  components  of  the  injector  flow: 

1.  A  distinct  head  end  region,  0  <  z  <  2 Rt,  where  the  bulk  flow  velocity  is  zero,  and  the  free 
surface  radius  is  equal  to  rv  as  shown  in  Fig.  1. 
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2.  A  nonzero  bulk  How  velocity  in  the  vortex  chamber,  in  the  region  2 Rt  <  z  <  Lv. 

3.  A  distinct  conical  convergence  section  connecting  the  vortex  chamber  to  the  nozzle,  spanning 
the  region  Lv  <  z  <  Lv  +  Lc. 

Item  1  tells  us  that,  now,  we  can  have  a  purely  standing  wave  only  in  the  head  end  region  of  the 
How.  Item  2  will  result  in  the  fact  that,  because  there  is  now  a  finite  bulk  fiow  velocity  in  the  vortex 
chamber,  the  lengths  of  the  waves  propagating  in  stream  wise  and  counter  stream  wise  directions 
will  differ  from  each  other,  which  leads  to  the  phenomenon  of  partial  standing  waves,  which  is  well 
described  in  Dean  and  Dalrymple  [17].  Lastly,  item  3  invokes  the  notion  of  a  smooth  variation  of 
bulk  fiow  velocity  and  boundaries  in  the  conical  convergence  region  up  to  the  point  where  the  free 
surface  radius  reaches  the  value  of  rn  in  the  nozzle.  The  nonuniform  fiow  causes  the  disturbance 
waves  to  both  refract  and  reflect  as  they  travel  through  the  fiow  transition. 

This  compound  problem  of  wave  refraction  and  reflection  is  difficult  to  attack  at  once.  However, 
there  is  a  simplifying  way  to  deal  with  it  by  saying  that  we  can  discretize  the  entire  transition  region 
into  short  cylindrical  sections,  in  wffiich  the  radii  of  solid  and  free  surface  boundaries  do  not  change, 
thereby  effectively  eliminating  refraction.  This  technique  is  very  similar  to  that  used  in  gravity 
wraves,  wffiere  the  classic  example  is  the  paper  by  O’Hare  and  Davies  [18].  But,  in  each  of  these 
short  sections,  wTe  need  to  knowr  the  local  w’ave  number.  It  can  be  showm  (see  [3,  Sec.  4.2])  that  the 
same  long  wTave  Eq.  (15)  can  be  used  in  this  case,  wfith  the  difference  that  the  local  values  of  uz , 
R ,  and  S  need  to  be  used.  In  turn,  the  local  bulk  fiowr  velocity  follows  simply  from  the  steady  state 
continuity  in  each  of  these  sections. 

We  will  start  the  analysis  first  by  considering  that  there  is  just  one  cylindrical  section  connecting 
the  vortex  chamber  to  the  nozzle.  This  wall  serve  as  a  good  platform  to  showr  the  main  features  of 
the  problem.  Then,  wTe  will  generalize  the  equations  used  in  this  simple  problem  for  further  using 
in  geometries  wffiere  more  transition  sections  are  considered.  The  analysis  will  be  concluded  writh 
an  algorithm  that  produces  solutions  for  such  general  geometries. 
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A.  CCRM  with  One  Cylindrical  Section  in  Transition 


For  now,  let  us  represent  the  transition  region  by  just  one  cylindrical  section  as  shown  in  Fig.  5. 
The  flow  parameters  corresponding  to  this  cylindrical  section  will  be  denoted  by  subscript  2  to 
indicate  that  it  is  located  next  to  the  nozzle  section.  Further,  let  the  radii  of  this  cylindrical  section 
be  given  as  the  average  between  the  corresponding  solid  boundary  and  free  surface  radii  of  the 
vortex  chamber  and  the  nozzle:  R2  =  (Rv  +  i?n)/2  and  r2  =  ( rv  +  rn)j 2.  We  have  four  sections  of 
the  How  connected  at  their  respective  discontinuity  boundaries:  the  head  end  region,  the  uniform 
vortex  chamber  region,  the  cylindrical  section,  and  the  nozzle.  Note  that  we  consider  the  nozzle  up 
to  the  point  where  we  assume  the  transition  ends,  which,  based  on  the  results  of  [3,  Sec.  3.4],  is 
located  at  z  =  Lv  +  Lc  0.5i?n,  see  Fig.  5. 


L  L  0.5R 


Fig.  5  Schematic  of  wave  reflection  and  transmission  for  CCRM  with  one  cylindrical  section 

Similarly  to  the  previous  model  ACRM,  if  an  incident  wave  M  exp  [i  ( ksz  —  cut)]  has  originated 
say  at  the  head,  it  will  generate  a  series  of  reflected  and  transmitted  waves  in  all  of  these  four 
flow  sections.  Fig.  5  shows  their  respective  directions  of  propagation  and  wave  numbers.  In  each 
section,  we  have  two  waves  traveling  in  opposite  directions.  We  can  calculate  their  wave  numbers, 
but  their  amplitudes  are  unknown.  Let  us  say  that  we  know  the  amplitude  of  the  original  incident 
wave,  M.  Since  there  is  a  purely  standing  wTave  at  the  head  end,  then  we  also  know  that  N  =  M 
and  /r7  =  —k8  (see  ACRM  for  more  discussion).  This  leaves  us  with  six  unknown  amplitudes,  A 
through  G.  There  are  three  matching  discontinuity  boundaries  that  can  relate  them  together,  with 
locations  at:  z  =  2Rt,  z  =  Lvy  and  2  =  Lv  +  Lc.  In  contrast  with  the  previous  ACRM,  at  each  of 
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these  boundaries,  in  addition  to  the  matching  of  the  fluctuating  mass  flow  rate,  we  will  also  assure 
the  matching  of  the  fluctuating  momentum.  This  will  accordingly  provide  six  equations  to  solve  for 
the  unknown  wave  amplitudes. 

Let  us  start  with  the  matching  of  the  fluctuating  mass  flow  rate  given  by  Eq.  (6)  at  z  =  2 i?*, 
z  =  Ly,  and  z  =  Lv  +  Lc.  Then  we  can  write  the  following: 


-Fel 


-Geu 


i  M 


.. ,  _  JkaL-  -  T2y-Eeik',L"  +  rv^Feik*L"  +  rv^-Geik*L'-  =  0 
k  3  k&  k6 


(16) 

(17) 


_  +  r  "  Dei>"(L„+Lc)  +  r  «  f.ei*4(t.+L.)  =  o  (18) 

k  i  k  2  ks  k  4 

Before  considering  matching  of  the  fluctuating  momentum,  let  us  introduce  the  following  coefficients  : 
Ka  =  -2T„u2zn  +  4uz„t„  — —  +  C2T  —  4 R2„rn  +  3 i?„r2) 

*1  Tn 

Kb  =  -2rnu2zn  +  4 +  C2-L  (*J  -  4 i?2r„  +  3i?„r2) 


Kd  =  —2tou2z2  +  4ua2r2  (i??  —  4i?2r2  +  3R-2T%) 

I<E  =  -2 7*2^2  H-  4^2 ^2  - —  +  C2-^  (^2  —  4i?2^2  +  3^2^) 

/14  r2 

A>  =  — 2r„52„  +  45,„r„  ~  ^  +  C21  -  4i?2r„  +  3fl„r2) 

*5  rv 

Kg  =  -2r„S2„  +  4a,.r.*6S?  ~  ^  +  C2i  (Jlj  -  4/?2r„  +  3i?„r2) 

*6  ry 

Then,  based  on  Eq.  (13),  the  matching  equations  for  the  fluctuating  momentum  at  z  =  2 i?*,  z  =  Lv, 
and  2  =  Ly  +  Lc  are  given  by: 


KFFeik*2R‘  +  KGGeik‘2R‘  =  C2  —  (R3V  -  4 i?2r„  +  3i?„r2)  (e-ik‘2R‘  +  eik‘2R‘)  M  (19) 

KDDeiksL“  +  KEEeik*L"  -  KEFe ik°L“  -  KGGeik‘L ”  =  0  (20) 

KAAeikAL~+LA  +  KBBeikAL»+LA  -  KDDeiks{L'’+LA  -  KEEeikAL-+LA  =  0  (21) 


We  can  rewrite  Eqs.  (1 6)— (1 8)  and  (19) — (21 )  in  a  matrix  form  as  shown  in  Fig.  6.  Solution  of 
this  matrix  equation  then  gives  the  dependence  of  the  wave  amplitudes  A  through  G  on  the  wave 
amplitude  of  the  original  incident  wave,  M. 
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What  information  can  we  extract  from  this  solution  with  regard  to  the  injector  resonance?  As 
before,  all  wave  numbers  in  the  matrix  equation,  Fig.  6,  depend  on  the  disturbance  frequency,  u. 
Therefore,  all  amplitudes  A  through  G  in  the  solution  depend  on  uj  as  well.  At  injector  resonance, 
because  the  film  thickness  fluctuation  in  the  injector  nozzle  should  be  maximum,  we  should  see  the 
greatest  magnitude  of  the  disturbance  wave  transmitted  into  the  nozzle,  which  is  represented  by 
the  amplitude  B ,  and  vice  versa,  we  should  see  the  smallest  magnitude  of  the  disturbance  wave 
traveling  back  into  the  vortex  chamber,  which  is  represented  by  the  amplitude  A.  Accordingly,  the 
frequency  at  which  \B\  is  maximum,  indicates  a  resonant  frequency. 

B.  Generalization  of  CCRM  to  Multiple  Cylindrical  Sections  in  Transition 

We  are  interested  in  having  multiple  cylindrical  sections  in  the  transition  region  in  order  to 
better  represent  a  conical  convergence  surface.  This  will  generate  a  larger  matrix  than  that  shown 
in  Fig.  6  for  finding  all  wave  amplitudes  involved  in  the  solution.  Let  the  solution  in  this  case  be 
represented  by  a  matrix  equation  X  •  a  =  Yy  where  X  and  Y  are  the  generic  matrices  similar  to 
the  left  hand  side  (LHS)  and  the  right  hand  side  (RHS)  matrices  in  Fig.  6,  and  a  are  the  generic 
wave  amplitudes.  To  illustrate  how  to  construct  the  matrices  X  and  Y ,  let  us  consider  again  the 
flow  setup  that  has  just  one  cylindrical  section  in  conical  transition  as  a  reference,  Fig.  7.  This  will 
serve  as  a  starting  platform,  from  where  the  solution  may  be  extended  further  on  to  more  transition 
sections. 


Head 


Fig.  7  Schematic  of  wave  reflection  and  transmission  for  generic  CCRM 
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Let  us  represent  each  flow  section  by  an  index  q ,  and  their  total  number  by  N.  Let  the  head 
end  section  be  denoted  by  the  index  q  =  N  +  1.  In  this  particular  case,  N  =  3  and  q  =  1  . .  .4.  Let 
the  first  index  under  the  wave  number  correspond  to  the  how  section  number,  and  the  second  index 
-  to  the  direction  of  wave  propagation:  1  -  for  upstream  traveling  waves,  and  2  -  for  dowmstream 
traveling  waves.  Moreover,  let  us  define  the  radial  boundaries  of  each  section  by  rq  and  Rqy  and  the 
discontinuity  boundaries,  at  which  the  how  sections  connect  to  each  other,  by  zq.  Finally,  note  that 
the  wave  amplitudes  a  are  counted  simply  in  a  sequence  1,  2,  3,  ... . 


Table  1  Values  of  rq,  Rq,  and  zq  in  generic  CCRM  with  one  cylindrical  section 


ri  =  rn 

Rl  =  Rn 

Z\  —  Lv  +  Lc 

T2  =  ( tv  +  r„)/2 

Rl  —  (Rv  +  i?n)/‘- 

Z2  =  Lv 

T  3  =  Tv 

R3  =  Rr 

z3  =  2  Ri 

r  4  =  rv 

a? 

II 

0? 

Z\  —  0 

By  comparing  Fig.  7  to  Fig.  5,  it  is  easy  to  determine  rq ,  Rq ,  and  zq ,  they  are  summarized  in 
Table  1.  Based  on  rq,  Rq ,  we  can  determine  the  bulk  how  velocity,  uzq ,  the  wave  numbers,  kq ,  and 
coefficients  Kq  (that  are  needed  for  the  momentum  balance,  see  page  18)  for  each  section  of  the 
how.  The  bulk  how  velocity  follows  from  the  one-dimensional  continuity  as 


,  __  RI-t2v 

*zq  —  Uzv  r 2  _  r2 
9  9 


(22) 


The  wave  numbers  for  both  upstream  and  dowTnstream  traveling  waves  can  be  found  based  on 
Eq.  (15)  as 


kq,l  ~  U- 


p-2  rq 

2  T* 


(23) 


kqt2  —  1 


X'C 

rI  -  r» 

U2  - 1 

zq  2  Ti 


(24) 


And,  by  using  Eqs.  (22)-(24),  we  can  determine  the  coefficients  K  as 

KqA  =  -2t„u2zv  +  4 uzqr„  +  C2 1  ( R*  -  4R2qrq  +  3Rqr2)  (25) 


kq,l 
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(26) 


Ki,2  =  -2r,<  +  “  +  C=i  _  4*’r,  +  gtf^) 

Kq,  2  T'g 

Now  we  can  rewrite  the  matrix  equation  shown  in  Fig.  6  in  terms  rq ,  zq,  and  A"g  as  given  in 
Fig.  8.  Similarly  to  the  previous  matrix  solution,  the  resonance  will  follow  from  this  solution  at 
frequencies  where  |a2|  is  at  its  maximum. 

An  algorithm  has  been  created  [3,  Sec.  5.5]  to  construct  the  matrices  like  that  shown  in  Fig.  8. 
The  matrix  is  divided  into  halves  resulting  from  mass  How  and  momentum  conditions  respectively 
and  constructed  following  standard  linear  algebra  techniques.  Solutions  have  been  obtained  for  a 
general  number  of  steps  as  described  in  the  following  section.  As  a  validation  of  this  algorithm,  we 
confirmed  that  it  replicated  the  results  obtained  in  the  previous  section  when  only  a  single  step  was 
present. 


VI.  Nondimensionalization  and  Baseline  Injector 

In  this  study,  we  will  use  the  fluid  density  p *,  nozzle  radius,  i?*,  and  the  mean  tangential  inlet 
inflow  velocity,  W*n9  as  dimensions.  The  dimensional  values  are  denoted  by  superscript  *.  Hence, 
for  the  parameters  that  we  will  use  further  on,  we  have: 

kt  =  kW'  C*  =  CWinRn 

‘n  n 

Since  eventually,  in  Part  II  of  this  research,  we  will  validate  theoretical  results  in  this  study 
against  Ahn’s  [10]  experimental  results,  the  baseline  injector  will  have  the  same  characteristics  as 
the  injector  used  in  his  experimental  testing  [10,  App.  B],  with  the  sizes  outlined  in  Table  2. 

The  baseline  steady  state  pressure  drop  and  convergence  angle  are:  Apinj  =  40.3  psi,  a  =  45°. 
Applying  the  methodology  described  in  Bazarov  [1,  Sec.  3.1],  we  can  calculate  the  steady  core  radii 
and  velocities,  Table  3. 


VII.  Results 

In  this  section,  the  results  will  be  presented  for  each  of  the  above  derived  models  of  the  wave 
reflection  and  resonance.  Both  ACRM  and  CCRM  are  compared  against  the  fundamental  condition 
for  resonance  shown  in  Eq.  (4).  All  results  presented  here  are  valid  for  the  baseline  injector  described 
in  Tables  2  and  3. 
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Table  2  Baseline  injector  geometry 


Parameter 

Dimensional 

Nondimensional 

Rn 

0.250  in 

1.0 

Rin 

1.125  in 

4.5 

Rt 

0.125  in 

0.5 

Rv 

1.250  in 

5.0 

Lt 

0.450  in 

1.8 

Ln 

1.000  in 

4.0 

Lv 

5.000  in 

20.0 

Nin 

4 

Table  3  Baseline  steady  core  radii  and  flow  velocities 


Parameter 

Dimensional 

Nondimensional 

The 

0.1T94  in 

0.7177 

Tv 

0.1794  in 

0.7177 

Tn 

0.2019  in 

0.8077 

win 

3.7596  m/s 

1.0 

Uzv 

0.1535  m/s 

0.0408 

Uzn 

10.8135  m/s 

2.8762 

A.  ACRM  Results 

From  Eq.  (4),  we  can  immediately  determine  the  resonant  frequencies.  A  practical  issue  results 
from  the  fact  that  the  ACRM  presumes  a  radial  step,  while  the  actual  injector  incorporates  a  conical 
convergence.  For  the  purposes  of  comparison,  we  shall  assume  that  the  step  is  located  at  the  center 
of  the  convergence  section,  which  implies  that  we  replace  the  vortex  chamber  length  in  Eq.  (4) 
with  the  length  Lv  +  Lcj 2,  where  Lc  is  the  length  of  the  convergence  section.  We  shall  apply  this 
methodology  throughout  as  we  assess  the  behavior  of  the  ACRM. 
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For  the  baseline  injector,  we  have  summarized  them  in  Table  4  for  the  first  five  modes.  Notice 
that,  due  to  this  calculation,  the  first  resonant  frequency  is  205.7  Hz,  to  which  we  will  refer  to  in 
Part  II  of  this  research,  when  comparing  with  the  experimental  first  resonant  frequency. 

Table  4  Resonant  frequencies  due  to  ACRM  (first  5  modes) 


Mode,  n 

CJQ 

Wq  (rad /sec) 

fS  (Hz) 

1 

2.2 

1292.2 

205.7 

3 

6.5 

3876.7 

617.0 

5 

10.9 

6461.1 

1028.3 

7 

15.3 

9045.5 

1439.6 

9 

19.6 

11630.0 

1851.0 

We  set  the  amplitude  of  the  incident  wave  at  the  head  end  to  unity,  D  =  1.  Then,  based  on 
Eq.  (14),  we  can  plot  the  amplitude  (absolute  value)  of  the  outgoing  wave,  \A\  versus  the  disturbance 
frequency,  Fig.  9.  Notice  that  the  peaks  of  |A|  are  located  exactly  at  the  same  frequencies  as 
was  predicted  by  the  fundamental  resonance  condition  (see  Table  4).  As  was  assumed  in  ACRM, 
these  peaks  indicate  the  frequencies,  at  which  the  mass  fiow  rate  pulsation  in  the  nozzle  becomes 
maximum,  and  hence  are  the  resonant  peaks. 


dimensional  f,  Hz 


Fig.  9  Amplitude  of  the  outgoing  wave  vs.  disturbance  frequency  in  ACRM-2 
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Also,  this  calculation  shows  that  the  value  of  \A\  varies  from  zero  to  9.45.  The  fact  that  \A\ 
shrinks  to  zero  at  some  frequencies  tells  us  that,  at  these  frequencies,  we  can  potentially  expect 
damping  of  flow  pulsations  in  the  injector.  And  the  fact  that  |A|  may  reach  to  values  as  large  as 
9.45  compared  to  the  value  of  D  =  1  signifies  that  the  amplitudes  of  the  outgoing  waves  may  be 
much  larger  than  the  amplitudes  of  the  original  incident  waves  in  the  vortex  chamber.  This  situation 
is  natural,  because  we  know  that  the  wave  length  in  the  nozzle  is  shorter  than  in  the  vortex  chamber, 
and  thus  the  wave  amplitude  has  to  be  larger  to  conserve  energy. 

B.  CCRM  Results 

From  the  way  the  CCRM  is  derived,  we  can  raise  a  question  of  how  many  short  cylindrical 
sections  there  should  be  in  the  transition  region.  Or,  in  other  words,  what  does  “short’-  mean  in 
relation  to  the  axial  dimension  of  the  transition?  We  can  deal  with  this  problem  as  follows.  We  will 
investigate  several  different  setups,  with  a  different  number  of  cylindrical  sections  in  the  transition 
in  each  case,  in  terms  of  the  resonant  peaks  that  each  of  them  produces.  Note  that  increasing  the 
number  of  cylindrical  sections  leads  to  larger  solution  matrices,  which  increases  the  computational 
time.  If  there  is  a  large  variation  in  the  answer  from  case  to  case,  we  should  be  looking  for  a 
converged  solution.  If  the  answer  does  not  change  much,  then  we  can  choose  the  setup,  that  is  the 
fastest. 

Let  us  first  consider  eight  cylindrical  sections  in  the  transition,  Fig.  10.  This  corresponds  to  the 
length  of  each  of  these  sections  equal  to  0.5 Rn.  The  amplitude  of  the  original  incident  wave  at  the 
head  end,  is  set  to  1  as  before.  For  this  setup,  the  amplitude  of  the  outgoing  wave,  do  (see  Fig.  7) 
versus  the  disturbance  frequency  is  shown  in  Fig.  11.  Notice  that  the  resonant  peaks,  now  located 
at  118  Hz  and  470  Hz,  are  different  from  the  first  two  resonant  peaks,  205  Hz  and  617  Hz,  predicted 
by  ACRM  (see  Table  4  and  Fig.  9).  This  immediately  indicates  that  the  resonant  characteristics 
of  an  injector  with  a  distinct  conical  convergence  section  do  differ  from  those  of  an  injector  with  a 
90°  step  transition.  Also,  by  looking  at  the  magnitude  of  |a2|  at  different  frequencies  (that  reaches 
the  value  of  1263  at  the  second  peak)  in  relation  to  the  magnitude  of  the  original  incident  wave, 
which  is  1,  one  may  wonder  why  the  former  is  so  much  larger.  The  answer  follows  from  the  fact 
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that  the  amplitude  of  the  upstream  traveling  waves  rises  to  infinity  as  we  move  closer  to  the  point 
where  the  transition  ends  and  the  uniform  nozzle  region  starts  (see  [3,  Chap.  4]).  This  fact  causes 
the  large  difference  between  the  amplitudes  of  the  original  incident  waves  and  the  outgoing  waves. 
Since,  after  all,  this  analysis  is  a  first  order  small  disturbance  analysis,  the  wave  amplitudes  are  not 
as  much  important  as  the  frequencies  where  they  peak  out. 

Now,  let  us  move  on  to  the  larger  number  of  cylindrical  section  by  decreasing  their  width  down 
to  5%  of  the  previous  width  0.5i?n,  The  width  is  denoted  as  w  in  this  investigation.  Fig.  12  shows 
the  new  locations  of  the  resonant  peaks  and  Table  5  summarizes  the  respective  resonant  frequencies 
along  with  the  sizes  of  the  solution  matrices  in  each  of  the  cases  considered. 

Table  5  Sensitivity  of  resonant  peaks  in  CCRM  to  cylindrical  section  width 


Width,  w 

1st  peak.  Hz 

2nd  peak,  Hz 

Solution  matrix  size 

0.500f?n 

118 

470 

10x10 

0.2501^ 

119 

463 

19x19 

0.l00Rn 

121 

484 

46x46 

0.050 Rn 

110 

462 

91x91 

0.025 Rn 

120 

455 

181x181 

In  Table  5,  we  can  see  that  the  resonant  peaks  are  moving  within  3%  of  the  baseline  values  as  we 
vary  w.  The  solution  matrices,  however,  grow  roughly  two  times  bigger  each  time  we  decrease  the 
wfidth  to  the  half  of  the  previous.  This  means  that  we  can  choose  the  baseline  case,  with  w  =  0.5i?n, 
for  further  calculations,  because  it  requires  the  least  computational  time. 

VIII.  Conclusions  and  Discussion 

In  this  paper  we  have  presented  possible  methods  for  accounting  for  the  disturbance  wave 
reflections  and  the  resonance  caused  by  them  in  the  rocket  swirl  injector.  We  have  considered  two 
models  differing  in  their  levels  of  complexity  depending  on  whether  or  not  we  consider  the  presence 
of  the  nozzle,  and,  if  we  do,  then  how  exactly  do  we  treat  the  connection  of  the  nozzle  to  the  vortex 
chamber:  either  through  a  sudden  step  discontinuity,  or  through  the  conical  convergence  section. 
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(a>  Normal  view 


00  Zoomed  in  view 


Fig.  10  Example  of  eight  cylindrical  sections  in  transition  considered  for  wave  reflections  in 
CCRM 


2* 


Fig.  11  Outgoing  wave  amplitude  vs.  disturbance  frequency  for  CCRM  with  eight  cylindrical 
sections  in  transition  (peaks  at  118  and  470  Hz) 


Fig.  12  Sensitivity  of  outgoing  wave  amplitude  to  cylindrical  section  width  in  CCRM  (resonant 
peaks  are  summarized  in  Table  5) 
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The  results  following  from  the  first  abrupt  contraction  resonance  model  (ACRM)  have  shown 
the  same  resonant  frequency  predicted  by  a  simple  wave  analysis,  Eq.  (4),  which  confirms  the 
existence  of  resonant  frequencies  within  a  swirl  injector.  Because  both  of  them  take  into  account 
only  a  sudden  step  discontinuity,  we  can  potentially  expect  that  both  of  them  can  be  used  for  the 
resonance  analysis  of  the  injectors  having  90°  sudden  convergence,  or,  steep  convergence  angles,  in 
general. 

From  the  results  of  conical  convergence  resonance  model ,  CCRM,  we  have  seen  that  the  am¬ 
plitude  of  the  outgoing  wave  can  be  much  larger  than  the  amplitude  of  the  incident  wave.  This 
is  attributed  to  the  fact  that  the  amplitude  of  the  upstream  traveling  wave,  at  the  point  where 
the  nozzle  entrance  transition  region  ends  and  the  uniform  nozzle  region  begins,  grows  to  infinity, 
as  follows  from  [3,  Chap.  4].  Also,  we  have  seen  that  the  resonant  frequencies  following  from  the 
CCRM  are  different  from  those  in  ACRM.  This  clearly  shows  that,  when  the  injector  has  a  distinct 
conical  convergence  between  the  vortex  chamber  and  the  nozzle,  its  wave  reflection  and  resonance 
characteristics  are  different  from  the  injector  having  the  90°  sudden  step  convergence. 

The  question  is  now:  can  we  trust  these  models  to  predict  the  injector  resonance?  Our  approach 
to  answer  this  question  will  be  through  setting  up  a  computational  BEM  model  in  Part  II  of  this 
study,  wThich  closely  replicates  the  boundary  conditions  used  in  the  analytic  models  here,  and  going 
through  the  parametric  study,  in  which  we  can  investigate  the  influences  of  such  parameters  as 
a,  Lv,  Ln,  Rv ,  and  Win  on  the  injector  response.  The  comparison  of  the  frequencies  where  the 
response  is  maximized  with  the  resonant  frequencies  predicted  by  the  analytic  models  in  this  part 
of  the  study  will  provide  an  indication  of  how  adequate  they  are.  Following  this  logic,  in  Part  II, 
we  will  be  presenting  both  theoretical  and  computational  results  in  parallel.  Ultimately,  we  shall 
compare  the  theoretical  and  computational  results  to  the  experimental.  While  there  are  limited 
data  for  this  comparison,  we  do  address  this  issue  in  Part  II  as  well. 

To  conclude,  it  has  to  be  noted  that,  in  such  reflection-refraction  problems,  where  the  transition 
is  approximated  by  the  cylindrical  sections,  in  addition  to  the  regular  linear  reflected  and  transmitted 
waves,  we  would  expect  another  type  of  waves  to  form,  which  would  damp  out  as  they  propagate 
far  away  from  their  respective  step  discontinuities.  But  usually,  the  required  analysis  applies  a 


30 


variational  approach,  a  classic  gravity  wave  demonstration  of  which  is  provided  in  Miles  [19].  In  this 
study,  we  have  ignored  these  waves,  as  the  application  of  variational  analysis  to  the  swirling  How 
is  difficult.  But  their  inclusion  could  potentially  make  the  reflection/ resonance  analysis  presented 
here  to  be  more  precise. 
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Appendix  B  -  On  the  Dynamic  Response  of  Rocket  Swirl  Injectors  Part  II. 

Nonlinear  Dynamic  Response 
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On  the  Dynamic  Response  of  Rocket  Swirl  Injectors 
Part  II.  Nonlinear  Dynamic  Response 
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Purdue  University.  701  W.  Stadium  Ave.,  West  Lafayette,  IN  47907 


Nonlinear  boundary  element  calculations  are  used  to  compare  and  contrast  results  from 
the  linear  theories  described  in  Part  I  of  the  study.  Parametric  evaluations  are  con¬ 
ducted  to  assess  the  influence  of  vortex  chamber  dimensions,  convergence  angle,  nozzle 
length,  injector  flow  rate,  and  pulsation  magnitude  on  dynamic  response  characteris¬ 
tics.  Resonant  frequencies  are  compared  against  the  linear  theory.  Overall  magnitude 
of  frequency  response  is  characterized  for  a  wide  range  of  injector  designs. 


Nomenclature 


C 

f 


dsgrid 


Lc 

Ln 

Lv 

m„ 

Nin 

Rin 

q 


=  angular  momentum  constant  (C  =  uer) 
=  disturbance  frequency 
=  grid  spacing 

=  length  of  conical  convergence  section 
=  length  of  nozzle 
=  length  of  vortex  chamber 
=  nozzle  exit  mass  flow  rate 
=  number  of  tangential  inlets 
=  inflow  radius  ( Rin  =  Rv  —  Rt) 

=  normal  velocity 
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Rn 

Rt 

Rv 


T  he 


T  n 

rv 

Ue ,  Ur ,  uz 

win 

a 

Apini 

nini 

<i> 

CJ 

prime  ()' 
bar  () 
hat  () 
star  ()* 


radius  of  nozzle 
radius  of  tangential  inlet 
radius  of  vortex  chamber 

steady  free  surface  radius  at  head  end,  z  =  0 

steady  free  surface  radius  in  uniform  nozzle  region 

steady  free  surface  radius  in  uniform  vortex  chamber  region 

velocity  components  in  circumferential,  radial,  axial  directions 

tangential  inlet  inflow  velocity 

angle  of  solid  wall  convergence 

total  injector  pressure  drop 

free  surface  deflection  away  from  5 

total  injector  response 

velocity  potential 

angular  disturbance  frequency 

fluctuation  value  of  parameter 

steady  state  value  of  parameter 

amplitude  of  fluctuation  of  disturbed  parameter 

dimensional  value  of  parameter 


I.  Introduction 

In  Part  I  of  this  study,  we  developed  two  analytic/linear  models  for  assessing  the  magnitude 
of  the  injector  masstiow  pulsation  induced  by  a  sinusoidal  pressure  disturbance  on  the  air-core  of 
a  classic  swirl  (simplex)  injector.  The  Abrupt  Convergence  Resonance  Model  (ACRM)  describes 
the  magnitude  of  injector  masstiow  response  assuming  a  radial  step  contraction,  while  the  Conical 
Convergence  Resonance  Model  (CCRM)  treats  the  contraction  from  the  vortex  chamber  to  the 
nozzle  using  an  arbitrary  number  of  small  radial  steps.  In  this  second  part  of  the  study,  numerical 
calculations  are  performed  to  assess  nonlinear  dynamic  response  and  to  compare  results  with  ACRM 
and  CCRM  predictions.  The  vehicle  for  these  computations  is  a  boundary  element  method  (BEM) 
that  preserves  surface  shapes  with  high  accuracy. 
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While  the  community  is  pursuing  nonlinear  behavior  in  numerous  free  surface  problems,  there 
are  a  comparitively  small  number  of  results  available  for  the  swirl  injector.  Notable  exceptions 
to  this  observation  include  prior  works  from  our  own  group  [1,  2],  as  well  as  full  Navier  Stokes 
computations  [3].  The  prior  efforts  from  our  group  utilized  an  axisymmetric  BEM  that  modeled  the 
inflow  tangential  channels  as  an  axisymmetric  slot  of  equivalent  inlet  area.  In  the  present  study, 
this  condition  is  relaxed  and  a  more  natural  axisymmetric  inflow  condition  is  applied  at  the  head 
end  of  the  vortex  chamber.  Utilizing  this  inflow  condition  also  provides  a  direct  comparison  against 
analytic  linear  models  described  in  Part  I  of  the  study.  Parametric  studies  are  conducted  to  assess 
the  effect  of  vortex  chamber  dimensions,  contraction  angle,  nozzle  length,  and  pulsation  magnitude 
over  a  range  of  frequencies.  Recent  experimental  work  [4,  5]  provides  a  strong  motivation  for  the 
present  work  as  it  clearly  depicts  evidence  of  a  resonance  condition.  As  evidence,  we  include  the 
experimental  result  from  this  work  depicting  the  frequency  response  of  spray  cone  angle,  Fig.  1.  The 
localized  peak  at  the  frequency  of  about  221  Hz  was  theorized  to  be  a  signal  of  a  resonant  mode. 


Pulse  frequency  (Hz) 

Fig.  1  Experimental  data  due  to  Ahn  [4,  Fig.  6.15]  on  spray  cone  angle  fluctuation  at  distances 
of  0.7  and  2  nozzle  diameters  from  the  nozzle  exit,  peak  observed  at  221  Hz 
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For  this  reason,  the  geometry  used  in  this  experiment  served  as  a  baseline  for  use  in  the  present 
study.  This  geometry  was  provided  in  Tables  3  and  4  of  Part  I  of  this  study. 

Continuing  the  discussion  on  that  experimental  peak,  one  would  wonder  what  would  the  existing 
Bazarov’s  [6]  theory  give  in  its  vicinity.  Figure  2  shows  the  answer.  There  is  no  peak  around  221  Hz, 
and  the  closest  peak  in  this  area  is  at  361  Hz,  then  comes  the  second  peak  located  further  in  the 
higher  frequencies,  at  726  Hz.  Our  motivation  is  to  assess  the  differences  between  the  Bazarov’s 
theory,  the  linear  theory  developed  in  Part  I  of  this  study,  and  fully  nonlinear  computations.  The 
following  section  provides  a  brief  review  of  the  computational  tool,  followed  by  results  of  parametric 
studies  and  conclusions  from  the  study. 


Fig.  2  Response  of  baseline  injector  calculated  using  Bazarov’s  [6]  linear  theory,  first  peak 
observed  at  361  Hz,  second  peak  observed  at  726  Hz 


II.  Model  Development 

Reference  [7]  provides  a  complete  description  of  the  basic  model  elements;  only  highlights  will 
be  presented  here  in  the  interest  of  brevity.  An  inviscid,  incompressible,  axisymmetric  how  is 
presumed  such  that  the  how  dynamics  are  governed  by  Laplace’s  equation,  V2<£  =  0  .  The  boundary 
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element  method  utilizes  an  integral  representation  of  this  equation  to  provide  a  connection  between 
velocity  potential  <j)  values  on  the  boundary,  the  local  geometry,  and  the  local  velocity  normal  to 
the  boundary,  q  =  d(j)/dny  as  follows: 

Ci<t>(fi)-\- j  qG  dT  =  0  (1) 

where  is  the  value  of  the  potential  at  a  point  r  j  ,  T  is  the  boundary  of  the  domain,  c,  is  the 

singular  contribution  when  the  integral  path  passes  over  the  “base  point, r  and  G  is  the  free  space 
Green’s  function  corresponding  to  Laplace’s  equation.  For  an  axisymmetric  domain,  the  free  space 
Green’s  function  can  be  expressed  in  terms  of  elliptic  integrals  of  the  first  and  second  kinds  and  is  a 
function  solely  of  the  instantaneous  surface  geometry.  For  this  reason,  a  discrete  representation  of 
Eq.  (1)  can  be  cast  as  a  linear  system  of  equations  relating  local  (j)  and  q  values.  In  the  discretization, 
both  (f)  and  q  are  assumed  to  vary  linearly  along  each  element,  thereby  providing  formal  second- 
order  accuracy  for  the  method.  Since  the  resulting  integrals  do  not  have  exact  solutions  in  this  case, 
Gaussian  quadrature  is  used  to  maintain  high  accuracy  of  integration  and  preserve  second-order 
accuracy  overall. 

While  this  governing  equation  is  linear,  nonlinearities  in  these  free  surface  problems  enter 
through  the  boundary  condition  at  the  interface.  With  regard  to  the  inflow  boundary  and  the 
solid  wall  boundary  we  can  set  the  normal  velocities,  q ,  exactly.  At  the  inflow  boundary,  they  may 
be  set  to  their  prescribed  values  as  a  function  of  time,  and  on  the  solid  boundary,  they  must  be 
zero  at  all  times.  Accordingly,  the  velocity  potentials,  0,  are  the  unknowns  on  these  boundaries. 
The  unsteady  Bernoulli  equation  provides  a  connection  between  the  local  velocity  potential  and 
the  surface  shape  at  any  instant  in  time.  Prior  formulations  [7]  have  provided  a  derivation  of  this 
result  suitable  for  implementation  in  a  Lagrangian  surface  tracking  environment.  For  the  swirling 
How,  modifications  are  required  to  account  for  the  centrifugal  pressure  gradient  created  by  the  swirl. 
In  the  present  study,  we  have  made  substantive  updates  to  both  inflow  and  free  surface  boundary 
conditions  and  for  this  reason  will  provide  additional  detail  here. 
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A.  Boundary  Conditions  and  Computational  Mesh 


In  prior  swirl  injector  simulations  by  Park  [1]  and  Richardson  [2],  the  tangential  inflow  was 
modeled  as  a  cylindrical  slot,  where  the  How  enters  through  its  upper  cylindrical  surface  and  proceeds 
radially  towards  the  vortex  chamber.  The  prior  inflow  condition  is  contrasted  with  the  current 
boundary  condition  in  Fig.  3.  We  choose  to  introduce  fluid  axially  at  the  head-end  of  the  vortex 
chamber  to  better  replicate  the  assumptions  of  the  theoretical  models  in  Part  I  of  the  study  [8], 
and  because  the  modeling  of  the  tangential  channels  is  inherently  three-dimensional.  Physically,  the 
inflow  plane  would  reside  just  downstream  of  the  tangential  inlet  holes  (a  distance  2 Rt  as  shown  in 
Fig.  3)  as  Bazarov  and  others  have  shown  that  there  are  slight  changes  in  the  free  surface  in  this 
region  due  to  the  conversion  from  a  pure  swirling  flow  to  one  that  also  has  an  axial  velocity.  The 


qin(t) 


Fig.  3  Schematic  of  current  BEM  flow  setup  compared  to  prior  models  by  Park  [l]  and 
Richardson  [2] 
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axial  velocity,  qiny  is  defined: 


Qin 


NinRjWl 

-  (n  (*))2 


where  Nin  is  the  number  of  tangential  inlets,  W in  =  1  is  the  steady  inflow  velocity  (per  the  nondi- 
mensionalization  we  employ),  r i  is  the  radius  to  the  free  surface  at  the  inlet  plane,  and  Rv  is  the 
radius  of  the  overall  vortex  chamber.  Since  the  free  surface  radius  r:  adjusts  to  balance  the  overall 
imposed  pressure  drop,  as  we  initialize  a  steady  flow,  qin  will  adjust  until  the  free  surface  is  no 
longer  undulating.  The  inflow  velocity  at  this  steady-state  condition  we  define  as  g,n.  Using  this 
average  value  over  the  entire  inflow  plane,  we  then  impose  a  sinusoidal  oscillation: 


Qin  (t)  =  Qin  +  Q',n  (*)  =  Qin  (1  +  <?o»c)sin  (wt) 


where  qosc  is  the  amplitude  of  oscillation. 

With  regard  to  the  gaseous  pressure  in  the  core,  linear  analyses  [9,  Sec.  6.3]  show  that  its  action 
may  be  noticeable  only  if  the  gas-to-liquid  density  ratio  is  on  the  order  of  10%.  Lastly,  we  will 
presume  that  the  overall  effect  of  the  surface  tension  presence  is  quite  weak  with  regard  to  the  wave 
dynamics  in  the  vortex  chamber  and  nozzle.  Let  us  generalize  these  conclusions  and  assume  that 
we  can  compltely  ignore  the  inertial  and  capillary  effects  of  the  core  further  on  in  this  study,  which 
allows  us  to  set  the  gas  pressure  in  the  core  to  zero,  pg  =  0. 

The  free  surface  boundary  condition  is  the  most  challenging  due  to  the  unsteady  and  swirling 
nature  of  the  how.  Physically,  the  free  surface  boundary  condition  can  be  thought  of  as  a  local 
pressure  balance  where  dynamic  and  static  pressure  is  interchanged  (potentially  with  capillary  pres¬ 
sures  as  well)  via  the  description  of  an  unsteady  Bernoulli  equation.  We  cannot  “see5'  the  swirling 
component  of  the  how,  which  is  perpendicular  to  axial/radial  how  plane,  yet  this  portion  of  the  how 
imposes  hydrostatic-like  pressures  that  must  be  reflected  in  the  overall  dynamic  behavior. 

Prior  implementations  due  to  Yoon  [10],  Park  [1],  and  Richardson  [2],  introduce  a  total  potential 
held,  </>*,  that  is  comprised  of  two  parts:  the  first,  which  describes  the  movement  only  in  the  radial 
and  the  axial  directions  and  the  second,  which  describes  the  movement  only  in  the  circumferential 
direction.  The  linearity  of  the  governing  Laplace’s  equation  permits  a  direct  superposition  of  these 
potentials.  The  resulting  Bernoulli’s  equation  of  these  authors  has  included  the  Rossby  number  as 


a  dimensionless  representation  of  the  swirl  level.  However,  for  a  dynamic  situation,  the  swirl  level, 
and  hence  the  Rossby  number,  can  vary  in  time  and  an  improved  treatment  was  required  in  order 
to  reflect  this  fact.  For  this  reason,  we  employ  a  slightly  different  formulation  in  the  present  work. 

Since  angular  momentum  must  be  conserved  in  a  potential  How,  for  an  axisymmetric  situation 
we  have: 


D  ,  x  DC  n 
Dt  v  J  Dt 


(2) 


wThere  D /Dt  is  the  material  derivative,  C  =  rug  is  the  angular  momentum  constant  associated  with 
the  How.  Since  C  must  be  invariant  in  a  Lagrangian  sense,  the  swirl  level  can  change  locally  as 
waves  convect  through  the  vortex  chamber.  The  overall  velocity  field  is: 


V<t>  =  [ur,u9,u2]  = 


d(j)  1  d(j)  d(j) 

'd<(>  C  d(j>' 

_dr'  r  dd'  dz_ 

_dr’  t  ’  dz_ 

(3) 


where  we  have  used  the  fact  that  the  circumferential  How  component  is  given  by  the  potential  free 


Next,  the  Lagrangian  derivative  for  the  whole  How  field  may  be  written  as 

D6  86 

-  ■  +  V0  •  V<£ 


(4) 


Dt  dt 

In  this  equation,  we  know  that  the  Eulerian  time  derivative  is  given  by  the  usual  unsteady  Bernoulli’s 
equation,  which,  since  we  neglect  the  gaseous  core  presence,  can  be  written  as 


(5) 


where  A  represents  the  steady  state  terms.  Without  the  loss  of  generality,  we  can  set  ^4  =  0  by 
incorporating  the  steady  terms  into  (j).  With  this,  we  can  rearrange  Eq.  (5)  to 


dd)  1 


The  second  term  in  Eq.  (4)  can  be  obtained  from  Eq.  (3)  as 

2 


(d<b\  (d6\ '  C2 

W-V0=  (/)  +(/)  +  — 

\or  J  \oz  J  r- 

Combining  Eqs.  (4),  (6),  and  (7),  we  have 

1  1  1  I"  fd6\2  ( dd>\ 2  C2 

=  --v<£-v^+v<£-v^=-v<£-V0=-U-_J  +(alj  +t t 


D4> 

7n 


(6) 


(7) 


(3) 
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To  get  the  Lagrangian  derivative  for  the  points  moving  in  the  2D  BEM  plane,  which  say  are 
represented  by  the  velocity  potential  <j>BEMi  we  can  subtract  the  Lagrangian  derivative  of  the 
swirling  component  of  the  How  from  the  Lagrangian  derivative  of  the  total  field, 

D^bem  _  D<p  D  ( C6 ) 

Dt  ~  Dt  Dt  (  ' 


The  Lagrangian  derivative  of  the  swirling  component  of  the  flow  can  be  written  as 


D(ce) 

Dt 


DC  n  DQ 
— — — 6  H-  C  — — — 
Dt  Dt 


(10) 


In  this  expression  we  know  that  the  first  term  is  zero  due  to  Eq.  (2),  and  for  the  second  term  we 
have 


de 

c 

~dl  + 

ury  ,  uz 
r 

de  Ids  dd 

dr  ’  r  dO1  dz 


C 

r2 


Then  we  can  rewrite  Eq.  (10)  as 


D  ( Cd )  _  C 2 

Dt  r5” 


Substituting  Eqs.  (8)  and  (11)  into  Eq.  (9)  we  obtain 

D<t>BEM  _l\(d<b\2  (d<i>\2  cA  C2  _  1  IV&A2  (d<f>Y  C°- 

Dt  2  \drj  \dzj  '  r2  r2  2  \\dr  J  \dz  J  r2 


(11) 


(12) 


To  be  consistent  with  the  <f)-q  notation  for  the  BEM  parameters  which  we  have  been  using  so  far, 
we  can  rewrite  the  last  equation  with  omitting  the  “BEMr  index, 


D<t>  _  l  \  fd0\2  ( d(f)\ 2  _  C 2 

Dt  2  \dr  J  ^  \dz  J  r 2 


(13) 


wThich  is  the  ultimate  form  employed  within  the  simulation. 

Fig.  3  also  shows  a  notional  computational  mesh  that  is  employed.  Because  the  surface  in  the 
vortex  chamber  undulates  with  time,  the  grid  on  the  inflow  boundary  is  allowed  to  compress  or 
expand  to  fill  the  instantaneous  distance  between  the  vortex  chamber  wall  and  the  current  head-end 
free  surface  location.  Fixed  grid  spacing  is  employed  on  solid  walls.  On  the  free  surface,  a  dynamic 
grid  is  utilized.  Using  a  Lagrangian  tracking  of  the  interface  points,  positions  are  updated  at  each 
instant  in  time.  A  spline  fit  is  then  applied  to  the  entire  free  surface  and  nodes  are  redistributed 
so  as  to  maintain  a  uniform  grid  spacing,  dsgrid>  Though  not  shown  in  the  figure,  ring-shaped 
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structures  are  formed  due  to  capillary  instabilities  on  the  conical  sheet  downstream  of  the  nozzle 
exit.  A  pinch  criteria  is  employed  to  shed  these  structures  when  nodes  on  opposing  sides  of  the 
ligament  get  within  a  specified  fraction  of  the  overall  mesh  spacing  [11].  Pinching  events  lead  to 
dynamic  changes  in  the  size  of  the  computational  mesh.  While  the  pinch  criteria  are  very  important 
when  attempting  to  establish  the  size  of  the  pinched  structures,  the  emphasis  of  the  present  study 
of  the  dynamic  massfiow  production  at  the  nozzle  exit  makes  this  criteria  of  lesser  importance  as 
demonstrated  in  subsequent  grid  function  convergence  studies.  Local  surface  slope  and  curvature 
are  computed  to  Ath  order  accuracy,  and  the  surface  is  evolved  in  time  using  a  4</l-order  Runge- 
Kutta  time  integration.  More  details  on  the  free  surface  treatment  can  be  found  in  Park  [1]  and 
Heister  [12]. 


B.  Injector  Response  Computation 

In  dynamic  system  evaluation,  the  injector  response  or  admittance  is  often  an  important  con¬ 
sideration.  The  dimensionless  response,  II ,  describes  the  level  of  massfiow  pulsation  created  by  a 
given  pressure  perturbation: 


H-  -  =  71ln 

inj  Ap'  ■ 

rmj 


AP,nj 


Because  in  our  current  BEM  model  we  are  analyzing  the  fiow  starting  from  the  left  end  of  the 
uniform  vortex  chamber  region,  or  from  z  =  2Rt  of  the  actual  injector,  and  the  tangential  inlets 
are  out  of  our  consideration  (see  Fig.  3),  then  the  total  pressure  drop  of  the  BEM  injector  is  equal 
to  the  pressure  drop  through  the  liquid  body  at  the  inflow  boundary.  The  average  pressure  drop 
required  to  establish  the  swirling  inflow  can  be  expressed: 


Apinj  =  P-Z- 


1 


1 


(14) 


where  we  note  that  the  free  surface  radius  at  the  inflow  boundary,  7*x,  is  a  function  of  time  when 
pulsations  are  present,  and  its  value  is  averaged  over  time  to  give  the  mean  pressure  drop.  It  can 
be  shown  [9,  Sec.  7.6]  that  the  unsteady  portion  of  the  pressure  can  be  expressed: 


&p’i  ni  =  P^rA 


(15) 
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Once  a  quasi-periodic  behavior  is  obtained,  the  value  for  r[  can  be  computed  as  the  difference  in 
maximum  and  minimum  surface  heights  over  a  given  period  of  the  oscillation. 

In  principle,  the  nozzle  masstiow,  mn,  can  be  computed  as  the  integral  of  the  axial  momentum 
tiux  across  the  film  at  the  exit  of  the  nozzle.  However,  as  the  computational  mesh  does  not  traverse 
the  film  in  this  region,  one  needs  to  compute  axial  velocities  at  presecribed  interior  points  as  was  done 
in  Richardson  [2].  However,  this  process  extends  computation  times  significantly  and  the  analysis  at 
the  curved  film  boundary  at  the  nozzle  exit  is  not  consistent  with  the  assumptions  employed  in  the 
analysis  in  Part  I  of  this  study.  For  these  reasons,  an  alternative  approach  was  taken.  Rather  than 
use  the  exact  exit  plane  to  measure  flowrate  histories,  we  use  a  point  a  distance  of  one  half  nozzle 
length  (0.5Ln)  from  the  exit  plane,  as  computations  show  that  the  film  is  very  nearly  parallel  with 
the  injector  axis  and  nozzle  wall  at  this  locale.  While  there  is  a  short  time  lag  between  this  location 
and  the  actual  exit  plane,  the  phase  shift  associated  with  this  lag  is  small,  as  the  velocities  in  the 
nozzle  are  large  and  typical  nozzle  lengths  employed  are  also  quite  small.  Moreover,  this  location 
permits  direct  comparison  against  the  analytic  results  developed  in  Part  I  of  the  study. 

To  avoid  the  inclusion  of  interior  points,  the  axial  velocities  on  either  side  of  the  film  at  the 
0.5 Ln  location  are  computed  using  4th  order  centered  difference  formulas  [9,  Sec.  7.6],  and  a  linear 
axial  velocity  profile  across  the  film  is  assumed.  In  prior  work  [2],  the  velocity  profiles  across  the 
film  were  shown  to  be  quite  linear  -  in  fact  the  velocity  on  either  side  of  the  film  are  nearly  identical 
due  to  the  inviscid  nature  of  the  assumed  flow.  This  process  permits  a  streamlined  evaluation  of 
the  nozzle  exit  flow  that  is  consistent  with  assumptions  employed  in  the  analytic  models  from  Part 
I  of  this  study. 


III.  Grid  Convergence  Study 

Four  different  mesh  spacings,  0.06,  0.07,  0.08,  and  0.09  were  evaluated  using  the  baseline  ge¬ 
ometry  described  in  Part  I  of  the  study.  Both  steady  and  unsteady  characteristics  were  evaluated. 
A  dimensionless  timestep  of  0.0005  was  used  in  all  simulations.  The  dimensionless  period  for  the 
highest  frequency  disturbances  to  be  evaluated  is  about  0.5  which  implies  that  we  take  about  1000 
timesteps  per  period  for  these  highest  frequencies. 
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Steady  state  surface  shapes  have  shown  to  be  converging  for  the  mesh  size  of  0.06  [9,  Sec.  7.7]. 
As  the  unsteady  nature  of  the  solution  is  of  primary  value  here,  the  overall  frequency  response  of  the 
baseline  injector  was  assessed  to  evaluate  mesh  sensitivity.  For  evaluation  of  unsteady  characteristics, 
a  sinusoidal  velocity  perturbation  of  amplitude  30%  of  the  mean  inflow  velocity  was  used  in  all  cases. 
Figure  4  shows  the  result  of  this  process  indicating  an  insensitivity  in  response  levels  to  the  four 
mesh  sizes  evaluated.  The  finest  mesh  spacing  of  0.06  was  selected  as  it  was  still  amenable  for  use  in 
the  computational  environment  available  for  the  study.  A  typical  run  time  on  a  2.4  GHz  processor 
was  roughly  60  hours  to  evaluate  response  at  a  single  frequency. 

IV.  Results 

Parametric  studies  were  performed  to  assess  the  influence  of  pulsation  magnitude,  injection 
conditions,  and  injector  geometry  on  its  dynamic  response.  Comparisons  are  made  against  the 
analytic  models  developed  in  Part  I  of  this  study.  Specifically,  we  compare  predicted  resonant  tones 
from  the  ACRM,  Abrupt  Convergence  Resonance  Model,  as  well  as  the  tones  predicted  from  the 
CCRM,  Conical  Convergence  Resonance  Model.  The  baseline  injector  geometry  described  in  Part  I 
of  this  study  serves  as  the  point  of  departure  for  all  calculations,  and  all  results  assume  a  30% 
pulsation  magnitude,  except  where  noted  in  the  pulsation  magnitude  assessment  in  the  following 
section.  In  total,  over  1000  individual  cases  were  evaluated  for  this  part  of  the  study. 

A.  Assessment  of  Nonlinearity  and  Choice  of  Pulsation  Magnitude 

The  effect  of  the  magnitude  of  the  velocity  pulsation  was  addressed  in  parametric  fashion  on 
the  baseline  injector  configuration.  Velocity  pulsation  magnitudes,  gosc,  equal  to  5%,  10%,  30%, 
and  50%,  of  the  mean  injection  velocity  were  considered. 

Figure  5  shows  the  magnitude  of  the  injector  responses  at  various  pulsation  levels.  Overall, 
the  response  is  very  insensitive  to  pulsation  level,  indicating  that  nonlinear  effects  are  of  secondary 
importance.  Similar  behavior  was  noted  in  assessing  dynamic  response  of  plain  orifice  [13]  and 
gas/liquid  coaxial  [14]  injectors. 
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Fig. 


Table  1  Response  of  free  surface  radii  T\  and  r2  to  various  pulsation  magnitudes  (baseline 
injector,  /*  =  476.3  Hz) 
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5% 

0.7200 

0.0007 

0.0981% 

0.8121 

0.0013 

0.1579% 

10% 

0.7200 

0.0014 

0.1961% 

0.8121 

0.0025 

0.3135% 

30% 

0.7199 

0.0042 

0.5868% 

0.8120 

0.0076 

0.9360% 

50% 

0.7199 

0.0070 

0.9772% 

0.8119 

0.0127 

1.5588% 

Figure  6  and  Table  1  show  the  time  histories,  and  the  average  and  fluctuation  values  of  the  free 
surface  radii  r i  and  r2  (where  the  latter  is  the  free  surface  radius  at  the  nozzle  center).  Perturbations 
in  the  surface  shape  are  rather  small  even  in  the  case  of  high  amplitude  masstiow  pulsations  as  the 
swirl  level  and  axial  velocity  adjust  for  the  flowrate  variation  with  little  response  visible  on  the 
actual  free  surface.  From  Table  1,  we  can  conclude  that  free  surface  fluctuation  increases  roughly 
linearly  with  the  increase  of  pulsation  magnitude.  Figure  7,  which  shows  the  wave  shapes  developing 
in  the  vortex  chamber  and  the  nozzle  at  different  pulsation  magnitudes,  at  time  t  =  100.  Within 
the  bulk  of  the  vortex  chamber  a  sinusoidal  wave  shape  is  present,  but  complex  shapes  evolve  in 
the  transition  and  nozzle  regions.  For  the  5  and  10%  pulsation  levels,  the  waves  are  very  small 
and  assessing  statistics  relative  to  masstiow  pulsation  magnitude  can  be  numerically  challenging. 
In  Fig.  7(b),  we  can  see  that  at  30%  pulsation  we  get  a  much  more  distinctive  wave  shape  in  the 
nozzle  than  at  5%  and  10%.  For  this  reason,  we  chose  this  pulsation  level  for  use  in  the  remainder 
of  the  parametric  studies.  At  the  same  time,  by  looking  at  the  time  histories  of  head-end  (rj) 
and  mid-nozzle  (r2)  radii  in  Fig.  6,  we  can  conclude  that  the  pulsation  can  be  well  described  as 
linear  at  all  considered  magnitudes.  A  question  may  arise  now:  why  do  we  get  a  linear  sinusoidal 
fluctuation  of  the  free  surface  at  pulsation  magnitudes  as  strong  as  30%  or  50%?  These  numbers 
seem  large,  however,  bear  in  mind  that  the  bulk  how  velocity  in  the  vortex  chamber  is  quite  small. 
Thus,  the  resulting  bulk  how  pulsation  may  be  considered  as  quite  weak,  which  accordingly  results 
in  the  linear  fluctuation  of  the  free  surface. 
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(a,)  Fluctuation  of  free  surface  radius  at  inflow  boundary 


Fig.  6  Stabilized  fluctuation  of  unsteady  free  surface  radii  as  a  function  of  pulsation  magnitude, 
g0»,  (baseline  injector,  /'  =  476.3  Hz,  pulsation  started  at  time  t  =  50) 
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B.  Conical  Convergence  Angle  Variation 


Variation  of  the  convergence  angle,  a,  does  not  have  a  great  influence  on  the  steady-state  shape 
of  the  free  surface  [9,  Sec.  3.1].  However,  as  changes  in  a  change  the  overall  length  of  the  injector, 
we  do  expect  for  this  parameter  to  effect  the  dynamic  response  of  the  injector. 

Figure  8  depicts  the  magnitude  of  the  frequency  response  of  the  baseline  injector  with  various 
convergence  angles.  Strong  response  is  evidenced  at  primary  resonance  and  at  the  second  harmonic 
of  the  resonant  frequency.  A  third  harmonic  appears  weakly  for  the  60  and  90  degree  cases,  but  is 
absent  at  lower  convergence  angles.  In  general,  increasing  a  tends  to  shift  resonant  peaks  to  higher 
frequencies  as  the  overall  injector  length  is  shortened  as  a  increases.  The  heights  of  the  response 
peaks  is  not  greatly  affected  until  a  reaches  the  30  degree  case.  The  long  convergent  section  in 
this  case  leads  to  some  destructive  interference  in  wave  patterns  thereby  dropping  and  broadening 
peaks  in  the  response  curve.  This  factor  may  also  explain  the  lack  of  higher  order  resonances  when 
a  values  are  below  60  degrees. 

Table  2  provides  a  comparison  of  computed  resonant  frequencies  with  those  calculated  from 
the  analytic  models  described  in  Part  I  of  the  study.  Since  the  abrupt  contraction  models  ACRM 
assume  a  vertical  wall  at  the  contraction  plane,  we  utilize  an  effective  length  that  includes  1/2  of 
the  contraction  length  to  evaluate  frequencies  and  resonance  conditions.  In  contrast,  the  CCRM 
does  retain  a  as  a  parameter.  In  general,  the  ACRM  is  in  better  agreement  with  BEM  for  the 
primary  resonant  peak  (Peak  1  in  Table  2),  while  the  CCRM  does  a  better  job  in  matching  the 
higher  harmonic  (Peak  2). 

Figures  9(a)  and  Fig.  9(b)  depict  the  response  levels  predicted  in  the  analytic  models  near  the 
resonant  conditions.  In  Fig.  9(a)  the  ACRM  model  shows  a  broad  peak  with  a  sinusoidal  shape 
in  the  region  near  resonance  while  BEM  results  show  a  much  sharper  peak.  The  ACRM  model 
does  a  poor  job  in  replicating  the  second  harmonic  as  was  noted  in  Table  2.  Figure  9(b)  show  the 
magnitude  predicted  by  the  CCRM  model.  These  amplitudes  show  an  even  broader  character  than 
the  ACRM  model,  as  a  wider  range  of  frequencies  lead  to  substantial  amplitude  waves  under  the 
CCRM  assumptions.  While  this  model  does  react  to  the  a  variation  in  the  same  manner  as  the 
BEM  has  shown,  the  peaks  have  the  overall  tendency  to  shift  to  the  higher  frequencies,  and  the 
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intensity  of  the  outgoing  waves,  a2 ,  does  not  exhibit  any  visible  trend. 


From  Table  2,  we  can  conclude  that  the  ACRM  agrees  with  the  first  BEM  peak  relatively  well 
(within  8%  in  all  cases)  but  does  a  rather  poor  job  in  replicating  BEM  results  for  the  second  peak. 
On  the  other  hand,  there  is  an  excellent  matching  of  the  second  peak  frequencies  of  BEM  with 
the  CCRM,  but  rather  poor  job  of  matching  the  primary  resonance  with  this  model.  The  CCRM 
delivers  large  refraction  waves  at  higher  convergence  angles  and  the  wave  interactions  are  highly 
complex  while  the  ACRM  delivers  simpler  waveforms.  For  this  reason,  the  explanation  as  to  why 
one  model  matches  one  peak  and  one  does  a  better  job  matching  another  peak  is  highly  complex  and 
a  topic  for  further  work.  Lacking  a  more  detailed  understanding  and  the  computational  resources 
required  for  a  nonlinear  calculation,  we  can  recommend  to  use  ACRM  for  primary  resonance  and 
CCRM  for  secondary  resonance. 


Table  2  Resonant  peaks  for  a  variation  cases  (based  on  Figs.  8  and  9) 


Q 

Peaks,  Hz 

BEM 

ACRM 

CCRM 

30° 

Peak  1 

178.6 

192.8 

116 

Peak  2 

416.7 

578.5 

448 

45° 

Peak  1 

193.5 

205.7 

118 

Peak  2 

464.4 

617.0 

470 

60° 

Peak  1 

208.4 

213.9 

120 

Peak  2 

494.1 

641.6 

471 

90° 

Peak  1 

232.2 

226.2 

n/a 

Peak  2 

529.8 

678.7 

n/a 

In  conclusion,  we  shall  take  a  moment  to  recognize  that  the  analytical  ACRM  and  BEM  values 
of  the  first  peak  for  a  =  45°  in  Table  2,  205.7  Hz  and  193.5  Hz,  fall  near  the  experimental  peak  of 
221  Hz  in  Fig.  1,  from  which  we  started  the  whole  discussion  in  this  part  of  the  study.  This  supports 
the  conclusion  that  the  resonance  condition  described  in  Eq.  (4)  of  Part  I  of  this  study  has  some 
merit  in  describing  the  frequency  where  a  locally  high  injector  response  could  be  expected.  We  also 
used  Bazarov’s  technique  to  analyze  this  experimental  condition  and  it  gives  a  resonance  at  a  much 
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higher  frequency  of  361  Hz  assuming  a  viscous  correction  value  of  0.1.  Ideally,  we  should  have  had 
a  second  experimental  peak,  whose  analytical  and  computational  values  we  already  have  in  Table  2, 
to  be  completely  certain  in  this  conclusion. 


C.  Vortex  Chamber  Length  Variation  (90°  convergence  angle) 

A  series  of  simulations  were  carried  out  to  assess  the  effect  of  vortex  chamber  length,  Lv  = 
5,10,15,19,  on  the  response  characteristics  of  the  baseline  injector  geometry  utilizing  a  sharp  step 
convergence  from  the  vortex  chamber  to  the  nozzle.  Figure  10  and  Table  3  show  respectively  the 
computational  BEM  results,  and  the  frequencies  where  the  BEM  and  ACRM  response  curves  for 
each  Lv  peak  out.  The  CCRM  cannot  be  used  here  because  a  cylindrical  section  needed  to  fit  into 
the  90°  step  transition  would  be  of  length  zero,  which  physically  collapses  the  CCRM  to  ACRM. 
We  can  see  that  the  peaks  shift  to  lower  frequencies,  as  the  vortex  chamber  length  is  increased.  An 
analogy  here  may  be  drawn  with  the  string  musical  instruments.  Smaller  size  instruments,  like  the 
violin,  produce  higher  pitch  sounds  than  the  bigger  size  instruments,  like  the  guitar,  or  contrabass. 
As  in  the  above  subsection,  this  is  attributed  to  the  fact  that  the  longer  vortex  chamber  naturally 
selects  /generates  longer  standing  waves  that  are  the  result  of  the  lower  pulsation  frequency,  and  vice 
versa.  We  can  look  at  this  from  the  mathematical  point  of  view  as  well,  if  we  rewrite  the  equation 
for  the  resonant  modes  from  Part  I  of  the  study, 


ljq  —  n 


2  Lv 


2  r*  ’ 


n=  1,3,5, ... 


(16) 


with  which  the  ACRM  peaks  are  calculated  in  Table  3.  From  this  equation,  it  becomes  clear  that 
the  values  of  the  resonant  frequencies  will  decrease,  as  the  vortex  chamber  length  is  increased.  Note, 
however,  that  this  equation  does  not  provide  any  information  about  the  amplitude  of  the  oscillation 
when  the  injector  is  at  resonance. 

The  Lv  =  5  case  provides  a  distinctly  different  character  with  a  broader  peak  and  no  evidence 
of  a  second  harmonic.  In  this  case,  complex  wave  shapes  are  obtained  as  the  nozzle  is  no  longer 
a  negligible  length  compared  to  that  of  the  vortex  chamber.  It  is  more  difficult  to  sustain  a  pure 
standing  wave  with  the  short  chamber,  but  substantial  amplification  is  still  available  in  a  variety  of 
complex  waveforms. 
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Continuing  the  discussion  on  the  oscillation  amplitudes,  in  Fig.  10,  we  can  also  see  that  the 
amplitudes  of  the  peak  responses  grow  larger,  as  the  vortex  chamber  becomes  longer.  The  explana¬ 
tion  of  physics  here  may  be  cast  in  terms  of  the  spring-damper  oscillator,  whose  vibration  energy 
is  conserved  at  all  times  and  is  given  by  E  =  Q.5moj2A2,  where  m  is  the  mass  of  the  oscillating 
body,  lj  is  the  oscillation  frequency,  and  A  is  the  oscillation  amplitude  (see  discussion  in  Kinsler  [15, 
Sec.  1.7]).  Notice  that  the  mass  has  the  power  of  1,  whereas  the  frequency  has  the  power  of  2.  Then, 
we  can  write  the  following:  (a)  the  increase  in  Lv  causes  proportional  linear  increase  in  the  mass 
of  the  liquid  body  in  the  injector’s  vortex  chamber;  (b)  we  also  know  from  above,  that  the  increase 
in  Lv  decreases  the  peak  frequency;  (c)  in  all  Lv  cases,  we  excite  that  liquid  mass  with  the  same 
energy,  which  is  the  kinetic  energy  of  excitation  0.5 qfn  at  the  inflow  boundary;  (d)  but,  if  we  have 
a  linear  mass  increase  from  (a)  and  a  quadratic  frequency  decrease  from  (b),  then,  to  conserve  the 
energy  of  oscillation,  the  amplitude  A  should  grow  in  the  above  formula  for  the  energy. 

Now,  let  us  compare  the  frequencies  where  the  peaks  are  located  in  the  BEM  responses  and  their 
analytic  counterparts.  In  Table  3,  we  can  see  that  the  ACRM  matches  the  first  resonant  peak  very 
well  for  Lv  values  at  or  above  10.  As  mentioned  above,  in  shorter  vortex  chambers,  the  wave  pattern 
cannot  be  described  as  a  simple  standing  wave.  Also,  it  is  clear  that  the  second  resonant  peak  is  far 
from  agreement,  which  follows  the  conclusion  in  the  previous  subsection  that  the  ACRM’s  cannot 

Table  3  Resonant  peaks  for  Lv  variation  cases  (based  on  Fig.  10) 


Lv 

Peaks,  Hz 

BEM 

ACRM 

Peak  1 

622.1 

754.1 

5 

Peak  2 

no  data 

2262.3 

10 

Peak  1 

395.9 

411.3 

Peak  2 

no  peak 

1234.0 

15 

Peak  1 

287.5 

282.8 

Peak  2 

610.8 

848.4 

19 

Peak  1 

232.2 

226.2 

Peak  2 

529.8 

678.7 
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capture  it.  On  the  other  hand,  the  matching  of  the  first  resonant  peak  strengthens  the  conclusion 
we  have  made  in  the  prior  subsection,  where  we  said  that  the  ACRM  may  be  successfully  used  to 
calculate  the  first  resonant  peak  in  steep  angle  geometries. 

D.  Nozzle  Length  Variation  (45°  convergence  angle) 

A  series  of  cases  were  run  at  nozzle  lengths  of  2,  4,  6,  and  8  using  the  other  dimensions  and 
fiow  conditions  from  the  baseline  injector  geometry.  Surface  shapes  were  evaluated  for  steady  inflow 
conditions  and  the  overall  film  shapes  differed  only  slightly  in  the  nozzle  entrance  region  for  the 
lengths  investigated.  The  shortest  nozzle  ( Ln  =  2)  did  show  a  slightly  different  free  surface  shape  in 
this  region,  but  the  longer  nozzles  were  nearly  identical.  The  computed  injector  responses  are  shown 
in  Fig.  11  and  their  peaks  summarized  in  Table  4  reveal  that  the  response  is  overall  insensitive  to 
the  nozzle  length.  Notice  that  the  Ln  =  2  response  curve  looks  same  as  the  others.  Which  tells  us 
that,  even  if  the  free  surface  in  the  nozzle  entrance  transition  region  was  different  in  this  case,  it  did 
not  affect  the  wave  reflection/transmission  characteristics  much. 

Theoretically,  the  deviation  of  the  steady  free  surface  in  the  Ln  =  2  case  from  the  rest  should 
affect  us  in  terms  of  the  CCRM,  where  the  point  wThere  the  transition  ends  is  one  of  its  inputs. 
Nonetheless,  based  on  the  knowledge  of  response  constancy  at  various  nozzle  lengths,  which  have 
learned  from  BEM  simulations,  and  for  simplicity,  we  will  assume  that  it  ends  at  the  same  +0.5i?n 
location  with  the  same  radius  equal  to  rn  as  the  other  free  surfaces.  This  allows  us  to  use  the  same 
baseline  injector  results  for  the  analytical  peaks  in  all  Ln  cases. 

Let  us  now  compare  the  analytical  and  the  BEM  peaks  with  each  other.  In  Table  4,  we  can  see 
that  the  BEM  computations  give  nearly  the  same  result  for  all  nozzle  lengths.  Similarly,  ACRM 
and  CCRM  values  do  not  vary  wfith  nozzle  length.  As  with  the  convergence  angle  study,  the  ACRM 
result  agrees  well  with  the  primary  peak  computed  in  BEM,  while  the  CCRM  better  matches  the 
second  peak.  In  general,  the  nozzle  length  has  little  influence  on  results  because  of  the  small  amount 
of  time  the  fluid  spends  in  this  region,  i.e.  characteristic  frequencies  in  the  nozzle  fiowpath  are  much 
higher  than  resonant  modes  determined  from  chamber  characteristics. 
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Table  4  Resonant  peaks  for  Ln  variation  cases  (based  on  Fig.  11) 


Ln 

Peaks,  Hz 

BEM 

ACRM 

CCRM 

Peak  1 

193.5 

205.7 

118 

2 

Peak  2 

464.4 

617.0 

470 

4 

Peak  1 

193.5 

205.7 

118 

Peak  2 

464.4 

617.0 

470 

Peak  1 

193.5 

205.7 

118 

6 

Peak  2 

470.3 

617.0 

470 

Peak  1 

193.5 

205.7 

118 

8 

Peak  2 

470.3 

617.0 

470 

E.  Vortex  Chamber  Radius  Variation  (45°  convergence  angle) 

A  series  of  simulations  were  conducted  with  differing  vortex  chamber  radii  (Rv)  while  keeping 
all  other  baseline  injector  geometry  and  how  conditions  the  same.  The  channel  inflow  velocity  was 
preserved  equal  to  the  value  described  in  Part  I  of  the  study  {W*n  =  3.7596  m/s)  which  implies  that 
the  overall  pressure  drop  was  varying  as  we  change  vortex  chamber  radius. 

Figure  12  shows  the  respective  how  geometries.  Notice  that  decreasing  the  vortex  chamber 
radius  leads  to  decrease  in  the  core  radius,  and  reduction  in  the  angle  of  the  conical  sheet  exiting 
the  nozzle.  This  is  due  to  the  fact  that  the  angular  momentum,  or  the  swirl  strength,  respectively 
becomes  smaller  as  we  have  preserved  the  inflow  velocity  for  all  cases.  In  turn,  this  allows  the 
swirling  fluid  to  descend  to  to  a  lower  core  radius.  And,  because  of  the  greater  inertia  in  the  axial 
direction  relative  to  the  baseline  case,  the  fluid  is  able  to  discharge  further  (see  similar  results  in 
Park  [1,  Fig.  4.15]). 

We  are  more  interested  in  the  behavior  of  the  free  surface  in  the  transition  region.  In  Fig.  12(c), 
we  can  clearly  see  that  the  transition  in  the  smaller  Rv  cases  starts  more  upstream  and  ends  more 
downstream  of  the  baseline  case.  In  this  study,  we  consider  the  following  approximations:  for 
Rv  =  3.75  case,  the  transition  starts  at  —  1.0i?n  and  ends  at  +1.0i?n,  for  Rv  =  2.50  case,  the 
transition  starts  at  — 1.5i?n  and  ends  at  +1.5i?„.  Accordingly,  the  corresponding  corrections  are 
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made  in  the  inputs  for  CCRM. 

Figure  13  and  Table  5  summarize  the  frequency  response  characteristics  for  the  three  Rv  values 
considered.  In  general,  the  vortex  chamber  radius  has  a  tremendous  effect  on  the  overall  response 
levels  as  larger  chambers  store  larger  amounts  of  fluid  and  become  more  susceptible  to  resonance 
pulsations.  As  we  decrease  the  vortex  chamber  radius,  the  peaks  shifts  to  the  left,  or  to  the  lower 
frequencies.  In  the  Rv  =  3.75  case,  we  start  to  see  the  third  harmonic  mode.  In  the  Rv  =  2.50  case, 
the  third  mode  is  very  apparent,  and  the  fourth  starts  to  appear.  Higher  order  harmonics  appear 
to  be  more  pronounced  (relative  to  the  primary  peak)  with  the  smaller  vortex  chambers. 

Mathematically,  the  shifting  of  the  peaks  to  the  left  can  be  clearly  attributed  to  the  decrease 
in  the  swirl  strength,  if  we  look  at  Eq.  (16),  which  shows  a  proportional  dependence  between  the 
resonant  frequency,  u.'0,  and  the  angular  momentum  constant,  C .  From  the  physical  point  of  view,  we 
can  say  that,  as  we  decrease  Rvy  the  relative  increase  of  the  flow  momentum  in  the  axial  direction, 
wThich  we  mentioned  above,  leads  to  natural  elongation  of  the  disturbance  waves,  which  in  turn 
causes  the  oscillating  flow  system  in  the  injector  to  “choose**  the  lower  resonant  frequencies. 

Let  us  now  take  a  look  at  the  results  of  the  analytic  resonance  models  in  Table  5.  We  can  see 
that  both  the  ACRM  and  the  CCRM  have  captured  the  above  BEM  trends  as  we  decrease  Rv  - 
shifting  of  the  peaks  to  the  lower  frequencies.  In  terms  of  the  actual  values  of  the  peaks,  however,  we 
can  see  in  Table  5  that  the  ACRM  performed  better  and  located  the  first  two  peaks  in  Rv  =  2.50  and 
Rv  =  3.75  cases.  Smaller  vortex  chambers  create  more  one-dimensional  flows  and  this  may  explain 

Table  5  Resonant  peaks  for  Rv  variation  cases  (based  on  Fig.  13) 


Rv 

Peaks.  Hz 

BEM 

ACRM 

CCRM 

2.50 

Peak  1 

70.2 

73.4 

no  peak 

Peak  2 

220.5 

220.3 

162 

3.75 

Peak  1 

134.2 

132.6 

no  peak 

Peak  2 

350.5 

397.7 

274 

5.00 

Peak  1 

193.5 

205.7 

118 

Peak  2 

464.4 

617.0 

470 
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why  the  ACRM  begins  to  agree  better  with  the  BEM  results  as  Rv  is  diminished.  In  general,  the 
CCRM  compares  more  poorly  as  Rv  is  diminished.  This  may  be  related  to  the  wrong  choice  of  the 
locations  where  we  have  assumed  that  the  free  surface  starts  and  ends  its  transitioning,  or  to  the 
fact,  that  in  the  smaller  Rv  geometries,  the  equations  of  the  long  wave  fluctuations  of  the  mass  How 
rate  and  momentum  are  not  accurately  representing  the  How  disturbances. 


F.  Flow  Rate  Variation  (45°  convergence  angle) 

In  this  study,  we  vary  the  steady  inflow  velocity  Win,  which,  one  should  note,  automatically 
causes  changes  in  the  incoming  mass  How  rate,  through  min  =  NinnRjWiny  and  changes  in  the  swTirl 
intensity,  through  C  =  WinRjn.  Varying  the  flowrate  had  only  minor  effects  on  the  overall  shape 
of  the  free  surface  as  higher  How  is  accompanied  with  higher  levels  of  swirl  in  all  cases.  Figure  14 
presents  the  frequency  response  for  several  flowrates.  Resonant  frequencies  decrease  with  decreasing 
flowrate,  and  the  amplitudes  of  the  peaks  decline  with  flowrate  as  well. 


Table  6  Resonant  peaks  for  Win  variation  cases  (based  on  Fig.  14) 


Win 

Peaks,  Hz 

BEM 

ACRM 

CCRM 

0.25 

Peak  1 

12.1 

12.9 

7 

Peak  2 

27.9 

38.6 

29 

0.50 

Peak  1 

48.4 

51.4 

30 

Peak  2 

119.1 

1.54.2 

118 

0.75 

Peak  1 

108.8 

115.7 

67 

Peak  2 

267.9 

347.1 

265 

1.00 

Peak  1 

193.5 

205.7 

118 

Peak  2 

464.4 

617.7 

470 

The  comparison  of  the  peaks  in  Table  6  reveals  that  the  CCRM  has  an  accurate  estimation  of 
the  second  resonant  peaks.  Note  that  the  value  of  Rv  in  this  study  was  favorable  for  the  this  model 
to  apply.  The  ACRM  still  does  reasonably  well  in  capturing  the  first  resonant  peak  frequency. 
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G.  Comparison  of  Resonant  and  Non-Resonant  Wave  Shapes 


For  this  study,  we  choose  a  90°  injector  with  the  vortex  chamber  length  of  Lv  =  19,  and  consider 
the  free  surface  motion  in  the  thin  region  around  the  mean  free  surface  radius  in  the  vortex  chamber 
over  a  period  of  the  How  pulsation.  During  this  time  interval,  we  will  first  present  the  results  for  the 
first  resonant  frequency  of  232.2  Hz  (see  Table  3),  and  then  the  results  at  a  non-resonant  frequency, 
which,  in  this  investigation,  we  have  chosen  to  be  357.2  Hz. 

Figure  15(a)  shows  the  mode  shape  of  the  resonant  wave  pattern  revealing  a  behavior  similar  to 
a  1/4  wave  oscillator  with  minimal  motion  near  the  node  at  the  head  end  and  maximum  motion  near 
the  nozzle  inlet.  This  1/4  wave  resonator  shape  was  postulated  in  development  of  the  Fundamental 
Condition  for  Resonance  and  ACRM  model  in  Part  I  and  explains  their  relative  success  in  replicating 
the  primary  resonant  peaks  in  the  parametric  studies.  We  should  note  that  there  is  some  motion 
evident  at  the  inflow  boundary  of  the  vortex  chamber  in  the  simulation  results  and  this  would 
explain  any  differences  between  BEM  and  ACRM  results. 

Now,  let  us  take  a  look  at  the  non-resonant  mode  shapes  in  Fig.  15(b).  The  situation  is  now 
completely  different.  We  do  not  have  a  distinct  node  or  an  antinode  and  the  standing  wave  behavior 
is  replaced  by  traveling  waves.  The  node  that  is  apparent  traveling  back  and  forth  around  the  center 
of  the  vortex  chamber.  This  reinforces  our  contention  in  Part  I  of  the  study;  when  the  injector  is  not 
at  its  first  resonant  mode,  the  wave  pattern  in  its  vortex  chamber  cannot  be  described  as  a  standing 
wave.  It  is  also  interesting  to  note  that  the  wave  pulsations  are  actually  larger  in  magnitude  than 
in  the  resonant  case. 
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Fig.  8  BEM  injector  response  sensitivity  to  conical  convergence  angle  variation 
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Fig.  10  BEM  injector  response  sensitivii 
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Fig.  11  BEM  injector  response  sensitivity  to  nozzle  length  variation  (45°  convergence  angle) 
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Fig. 
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Fig.  14  BEM  injector  response  sensitivity  to  steady  tangential  inlet  inflow  \ 


(a.)  First  mode  resonant  wave  pattern  (/*  =  232.2  Hz) 


(b)  Non-resonant  wave  pattern  (/*  =  357.2  Hz) 


Fig.  15  Resonant  and  non-resonant  free  surface  wave  envelopes  in  swirl  injector’s  vortex 
chamber  (a  =  90°,  Lv  =  19,  t  =  67.6950 .  .  .  70.0950) 


V.  Conclusions  and  Discussion 


The  dynamics  of  the  classical  swirl  injector  are  highly  complex  in  that  the  free  surface  in  the 
vortex  chamber  can  sustain  complex  wave  shapes  due  to  a  combination  of  reflections  and  refractions 
from  the  convergent  surface  that  forms  the  nozzle.  Large  dynamic  response  can  be  generated 
leading  to  dimensionless  massfiow  pulsations  many  times  greater  than  dimensionless  pressure  drop 
perturbations.  A  resonant  condition  has  been  shown  to  exist  wherein  the  injector  behaves  as  a  1/4 
wave  oscillator  thereby  creating  large  amplitude  pulsations  at  the  nozzle  that  lead  to  corresponding 
pulsations  in  massfiow.  The  theoretical  prediction  for  this  resonant  condition  (16)  agrees  quite  well 
with  the  fully  nonlinear  calculations  over  a  wide  range  of  injector  designs  and  flow  conditions. 

We  have  started  this  part  of  the  study  with  the  notion  that,  for  the  baseline  injector,  we  have 
one  experimental  data  point  for  spray  cone  fluctuation  at  the  pulsation  frequency  of  221  Hz.  We 
also  know  from  Part  I  of  this  study  that  the  ACRM  predicts  the  first  resonant  frequency  of  205.6  Hz 
while  the  nonlinear  BEM  simulations  give  a  frequency  of  193.5  Hz.  On  the  other  hand,  Bazarov’s  [6] 
response  curve  does  not  show  any  extremum  in  the  area  of  those  frequencies.  Based  on  these  limited 
data,  there  is  some  evidence  that  the  resonance  condition  described  in  this  work  has  merit.  Clearly 
it  would  be  desirable  to  compare  the  model  against  other  datasets.  Bazarov’s  book  contains  some 
high  frequency  data  he  took  in  the  1970’s,  but  the  description  of  the  injector  geometries  for  those 
tests  is  incomplete.  Hopefully,  the  present  work  will  motivate  some  more  fundamental  experiments 
for  which  we  can  compare  the  techniques  discussed  in  this  study. 

The  vortex  chamber  radius  has  the  most  prominent  effect  on  injector  response  with  smaller  Rv 
values  leading  to  smaller  response  levels  at  lower  frequencies.  However,  resonances  appear  to  be 
more  pronounced  as  Rv  is  decreased.  The  vortex  chamber  length  also  has  a  strong  influence  on  the 
levels  of  injector  response  with  longer  chambers  showing  higher  response  levels  at  lower  frequencies. 
At  very  small  chamber  length  ( Lv  =  5  case)  we  saw  a  fundamentally  different  character  of  response 
with  sharp,  well  defined  peaks  being  replaced  with  a  broader  peak.  The  chamber  massfiow  rate 
also  showed  strong  influence  on  the  overall  response  with  higher  flowrates  leading  to  larger  response 
levels  at  higher  frequencies. 
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The  nozzle  convergence  angle  (a)  and  nozzle  length  (Ln)  had  a  much  smaller  influence  on  overall 
injector  response.  Decreasing  a  tended  to  decrease  the  magnitude  of  resonant  peaks  slightly  and 
shift  the  frequencies  to  lower  values  as  the  overall  tiowpath  length  increased  with  this  adjustment. 
At  low  values  of  a  the  peaks  became  less  sharp  as  more  complex  wave  patterns  were  generated  due 
to  the  long  convergent  section.  The  nozzle  length  had  almost  no  influence  on  the  results  as  for 
practical  designs  even  long  nozzles  are  short  compared  to  the  vortex  chamber  length  and  the  fluid 
residence  time  in  the  nozzle  is  short  due  to  the  high  velocities. 

The  comparison  of  computational  and  analytic  results  is  somewhat  mixed.  The  ACRM  seems 
to  do  a  resonable  job  in  predicting  the  primary  resonance  peak  for  most  geometries.  In  principle,  the 
simple  resonant  frequency  relation  (16)  provides  the  designer  with  a  simple  mechanism  to  predict 
this  primary  resonance  if  one  uses  the  length  to  the  mid-point  of  the  convergent  section  as  the 
equivalent  termination  of  the  vortex  chamber.  Predicting  the  second  harmonic  (second  peak  in  the 
response  curve)  is  more  problematic  as  neither  the  ACRM  or  CCRM  results  compared  consistently 
well  with  the  BEM  calculations  over  the  parameter  space  investigated.  The  theoretical  resonance 
equation  (16)  shows  that  the  second  harmonic  should  be  at  a  frequency  3  times  that  of  the  primary 
harmonic  (i.e.  n  =  3  in  Eq.  (16)).  If  one  regards  the  BEM  results  as  exact  and  looks  at  the  ratio  of 
the  ACRM  prediction  with  BEM  calculations,  the  actual  n  values  vary  over  the  range  of  2.2  <  n  <  3 
based  on  results  in  Tables  2-6.  As  the  vortex  chamber  becomes  small,  we  attain  n  values  closer  to 
the  theoretical  value  of  3.0. 

While  the  CCRM  utilizes  more  physics  in  that  it  includes  a  momentum  balance  as  well  as  a 
mass  balance,  it  does  not  replicate  the  primary  resonance  frequencies  well  and  had  inconsistent 
agreement  with  the  second  harmonics  computed  from  BEM.  We  view  this  as  an  area  of  further 
study  as  it  is  not  well  understood  why  such  a  situation  should  exist.  While  the  CCRM  does  seem 
to  have  some  success  in  replicating  the  second  resonance  peak  frequency  for  a  number  of  different 
chamber  geometries  studied,  there  were  also  cases  with  substantial  disagreement  that  cannot  be 
fully  explained.  In  general,  analytic  models  will  have  difficulty  when  wave  patterns  become  more 
complex  due  either  to  geometric  or  input  signal  variations  that  lead  to  additional  complexity. 
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Appendix  C  -  Experimental  Study  Swirl  Injector  Dynamic  Response  using  a  Hydro - 
Mechanical  Pulsator 

Benjamin  Ahn  and  Maksud  Ismailov  and  Stephen  D.  Heister 

The  dynamic  response  of  a  classical  (simplex  style)  swirl  injector  has  been  studied 
experimentally  using  a  super-scale  transparent  model  with  water  as  the  working  fluid.  A 
unique  mechanism  was  developed  for  imparting  controlled  perturbations  to  the  injector  inlet 
massflow  by  successively  blocking  and  opening  tangential  inlet  flow  passages  using  a  rotating 
cap  over  the  inlet  ports.  Two  vortex  chamber  designs  (long  and  short)  were  evaluated  to  assess 
the  effect  of  this  important  design  variable.  High  speed  imaging  of  the  spray  cone  and  air 
core/liquid  interface  inside  the  vortex  chamber  was  used  to  assess  dynamic  behavior  at 
frequencies  up  to  500  Hz.  Resonant  conditions  were  detected  in  both  designs  and  both 
measurements  gave  similar  frequencies  for  the  resonant  peak.  The  resonant  peak  was 
compared  against  recent  theory  due  to  Ismailov11  and  results  compare  well  only  when  the 
theory  is  adjusted  to  account  for  potential  water  hammer  effects  induced  by  the  rotating  cap. 


Nomenclature 

R ullage 

= 

ullage  pressure 

R 'manifold 

= 

manifold  pressure 

Dn 

= 

diameter  of  nozzle 

Rn 

= 

radius  of  nozzle 

Rin 

= 

radius  to  centerline  of  inlet  channel  (Fig.  1) 

Rt 

= 

radius  of  tangential  inlet 

Rvc 

= 

radius  of  vortex  chamber 

rvc 

= 

radius  of  air  core  in  vortex  chamber 

Rvc 

= 

vortex  chamber  length 

U 

= 

effective  vortex  chamber  length  including  lA  nozzle  contraction  length 

vt 

= 

velocity  in  tangential  inlets 

0 

= 

total  spray  angle  fluctuation  amplitude 

@avg 

= 

average  total  spray  angle 

eta 

= 

core  disturbance  amplitude  (half  air  core  diameter  amplitude  over  mean  diameter) 

I.  Introduction 


Swirl  injectors  have  a  wide  array  of  application  to  both  airbreathing  and  rocket  combustors  due  to  their  ease  of 
manufacture  and  excellent  spray  production  characteristics.  The  theory  for  steady  operation  of  these  devices  is  well 
established  (at  least  for  ideal  inviscid  fluids)  and  results  from  these  analyses  have  been  calibrated  and  correlated 
against  numerous  experimental  data1.  In  contrast,  the  dynamic  characteristics  of  these  devices  have  been  much  less 
studied2,3,4,5.  There  is  ample  motivation  to  study  the  dynamic  behavior  as  in  propulsion  and  combustion  applications, 
there  are  numerous  opportunities  for  the  injector  to  resonate  with  processes  in  the  combustor  or  downstream  devices 
leading  to  potentially  damaging  consequences. 

The  presence  of  the  vortex  chamber  (Fig.  1)  which  interacts  directly  with  the  downstream  environment  due  to  the 
presence  of  a  gas  or  air  core  at  its  center,  leads  to  unique  dynamic  properties  for  this  device.  The  strong  internal 
feedback  coupling  of  the  flow  rate  and  pressure  fluctuation  across  different  parts  of  the  swirl  injector  lead  to 
fluctuations  in  the  volumetric  liquid  flow  rate  at  the  exit  of  tangential  inlets  into  the  vortex  chamber,  which  is  delayed 
by  certain  phase  shift.  These  pulsations  are  manifested  in  waves  that  traverse  the  air  core  of  the  vortex  chamber,  grow 
in  amplitude  as  the  film  thins  in  the  nozzle,  and  exit  into  the  spray  cone/film  formed  in  the  near  exit  region.  The 
surface  waves  lead  to  fluctuations  in  the  nozzle  flow  which  will  in  turn  influence  film  thickness,  spray  formation,  and 
spray  angle.  In  addition,  other  properties,  such  as  circumferential  velocity  and  pressure  drop  fluctuations  in  the  liquid 
vortex,  will  all  vary  due  to  pressure  fluctuations  across  the  injector.  All  these  coupling  effects  will  result  in  unsteady 
drop  size  distribution  and  spray  angles. 

Most  of  the  prior  knowledge  in  injector  dynamics  has  been  contributed  by  Dr.  Vladimir  Bazarov,  whose  studies  at 
Moscow  Aviation  Institute  span  much  of  the  past  four  decades2,3,6,7.  More  recently,  investigations  in  South  Korea8 
and  the  United  States9'13  have  amplified  on  the  initial  body  of  work.  Bazarov’s  initial  theory  shows  some  opportunity 
for  injector  resonance  at  intermediate  frequencies,  but  recent  work  of  this  group  has  established  a  simple 
relationship9,11  to  approximate  resonant  frequencies  using  the  analogy  of  waves  entering  a  harbor  of  a  fixed  volume. 

Experimental  work  in  this  field  is  challenging  as  most  of  the  dynamic  character  of  interest  occurs  at  very  high 
frequencies  of  hundreds  if  not  thousands  of  Hz.  Hydromechanical  pulsators  tend  to  drop  off  in  pulsation  amplitude  at 
these  frequencies,  but  they  have  been  used  successfully  in  a  number  of  studies8,14.  Piezoelectric  devices  are  capable  of 
much  higher  frequencies,  but  pulsation  amplitudes  (stroke  lengths)  tend  to  be  somewhat  limited.  One  example  of  high 
amplitude  modulation  at  high  frequencies  is  due  to  Anderson15. 

In  the  present  study,  it  was  decided  to  use  super-sized  injector  in  order  to  reduce  the  frequencies  over  which 
response  might  be  detected.  In  an  attempt  to  maximize  the  flow  disturbance  due  to  pulsation,  a  design  was  developed 
to  physically  block  tangential  inlet  channels  using  a  rotating  cap  with  multiple  holes  around  its  periphery.  The 
objective  of  the  work  was  to  quantify  the  dynamic  response  over  a  range  of  frequencies  and  flow  rates  for  two  distinct 
injector  designs.  The  following  section  provides  a  description  of  the  test  apparatus  and  facility,  and  results  and 
conclusions  are  provided  subsequent  to  this  discussion. 

II.  Experimental  Apparatus  and  Methods 


Swirl  Injector  and  Pulsator 

Figure  1  provides  a  schematic  representation  of  the  super-sized  swirl  injector  design  utilized  in  the  study.  Two 
units,  using  short  and  long  transparent  vortex  chambers,  were  fabricated  to  use  a  common  nozzle  design  with  using  a 
45  degree  half-angle  converging  section  as  depicted  in  Fig.  1.  Table  1  provides  a  summary  of  major  dimensions  in  the 
two  test  articles. 

The  pulsator  consists  of  a  shaft  (a)  driven  by  an  electric  motor,  which  rotates  a  dynamic  inlet  cap  (e)  that  has  16 
tangential  holes  of  the  same  diameter  as  the  holes  of  the  static  inlet  cap  (d)  as  shown  in  Fig.  2.  As  the  dynamic  inlet 
cap  (e)  rotates  about  the  shaft,  the  holes  in  it  align  periodically  with  the  holes  in  the  static  inlet  cap  as  shown  on  the 
right  side  of  Fig.  2,  thus  allowing  the  liquid  to  flow  into  the  vortex  chamber.  When  two  sets  of  holes  meet  and  start  to 
overlap,  a  mutual  cross-sectional  area  starts  to  form  and  reaches  a  maximum  value  when  two  sets  of  holes  are  fully 
aligned.  The  area  then  starts  to  decrease  as  the  dynamic  inlet  cap  continues  to  rotate.  It  is  the  change  in  this  mutual 
cross-sectional  area  that  gives  the  pulsating  effect  in  this  design.  Different  numbers  of  holes  in  the  dynamic  inlet  cap, 
as  well  as  variable  rotational  speeds,  generate  pulsation  frequencies  ranging  from  0  Hz  to  500  Hz,  at  a  maximum 
manifold  pressure  of  80  psi. 

The  swirl  injector  and  the  pulsator  are  assembled  into  the  manifold  (g),  which  includes  a  pressure  tap  (k)  on  the 
manifold  sidewall  aligned  with  the  tangential  swirl  chamber  inlets  with  its  axis  lying  in  A-A  plane.  A  Druck  pressure 
transducer  with  a  response  frequency  of  2.5  kHz  is  used  to  measure  the  pressure.  Sampling  rate  for  pressure 
measurement  was  1000  samples/sec.  The  entire  vortex  chamber  is  made  from  transparent  cast  acrylic  to  observe  the 
air  core  and  film  thickness  inside  the  swirl  chamber,  and  the  pulsator  and  manifold  are  made  of  304  stainless  steel. 


Tangential  inlet  orifice 


Vortex  Chamber 


Table  1:  Swirl  Injector  Specifications 


Parameter 

Short  Injector 

Long  injector 

Radius  of  nozzle 

Rn 

(inch) 

0.250 

(mm) 

6.350 

(inch) 

0.250 

(mm) 

6.350 

Inflow  radius 

Rin 

1.125 

28.575 

1.125 

28.575 

Radius  of  tangential  inlet 

Rt 

0.125 

3.175 

0.125 

3.175 

Radius  of  vortex  chamber 

Rvc 

1.250 

28.575 

1.250 

28.575 

Length  of  tangential  inlet 

Lt 

0.450 

11.430 

0.450 

11.430 

Length  of  nozzle 

Rn 

1.000 

25.400 

1.000 

25.400 

Length  of  vortex  chamber 

Rye 

5.000 

127 

6.125 

155.575 

Dynamic  iti let  cap  (rotating  counter  clockwise) 
Injector  inlet  cap  (stationary) 

t 


Fully  aligned  Section  A- A 


To  electric  motor 

* - 


Figure  2:  Left  -  Schematic  drawing  of  swirl  injector  pulsator:  a)  shaft,  b)  vanes,  c)  swirl  chamber,  d) 

injector  inlet  cap,  e)  dynamic  inlet  cap,  f)  manifold  cap,  g)  manifold,  h)  bearings,  i)  water  supply  inlet,  j,k) 
pressure  transducer  taps.  Right  -  View  of  cross  section  A- A,  where  holes  of  the  dynamic  inlet  cap  and  injector 

inlet  cap  are  aligned. 


Test  Facility  Design 


The  swirl  injector  experiment  (Fig.  3)  was  housed  in  the  Two  Phase  Flow  laboratory,  Neil  Armstrong  Hall  of 
Engineering,  at  Purdue  University.  The  water  storage  tank,  which  stores  120  gallons  of  water,  is  certified  for  150  psia 
maximum  ullage  pressure.  A  1-1/2  inch  manual  ball  valve  is  used  as  the  run  valve.  The  air  which  is  supplied  by 
building’s  compressor  and  regulated  by  a  300  psig  manual  pressure  regulator  is  used  to  pressurize  the  water  tank. 
Water  stored  in  a  pressurized  air  tank  is  transported  to  the  test  article  via  a  1-1/2  inch  flex-hose,  which  is  followed  by  a 
0.5  inch  flex  hose.  Water  fills  the  manifold,  enters  through  the  pulsator,  and  exits  thorough  the  nozzle.  Test  articles, 
optical  instruments,  and  camera  are  mounted  on  a  10  by  4  foot  optics  table.  The  spray  discharging  into  ambient 
environment  is  collected  using  a  clear  15  by  15  inch  square  tube  before  entering  a  dump  tank  for  it  to  be  drained.  The 
dump  tank  is  a  30  gallon  plastic  tank  with  a  0.75  inch  drain  valve.  For  flow  visualization,  a  high  speed  camera  is 
placed  downstream  of  the  nozzle  exit,  opposite  to  a  concentrated  light  source.  A  6x24  inch  1  KW  Altman  Lightning 
Co  stage  light  is  used  as  a  back  light.  A  Fresnel  lens,  6  inch  in  diameter,  0.06  inch  thick,  is  used  to  bend  and  focus  the 
rays  to  form  a  single,  concentrated  beam  of  high  intensity  light  from  the  stage  light.  The  focused  light  then  passes 
through  a  0.5  inch  thick  acrylic  diffuser  plate  before  reaching  the  plane  of  the  test  article.  The  experiment  was  tested 
with  a  horizontal  flow  with  the  assumption  that  the  gravity  forces  were  unimportant  due  to  the  internal  flow 
acceleration  from  vortex  chamber  to  nozzle  being  much  greater  than  the  gravity.  Reference  10  provides  additional 
details  on  the  facility  design  and  operations. 


Figure  3:  Schematic  of  experimental  apparatus. 

Image  Acquisition  and  Processing 

Images  are  captured  using  a  high  speed  camera  and  a  commercial  software  video  package16.  A  36  mm  zoom  lens  is 
used  to  enhance  the  images,  by  varying  the  focal  length  throughout  the  experiment.  Prior  to  each  set  of  tests,  images 
are  calibrated  (pix/inch)  by  placing  a  25  lines-per-inch  resolution  plate  at  the  injector  axis.  The  high  speed  camera 
settings  for  spray  and  air  core  acquisition  are  tabulated  in  Table  2.  The  procedure  for  detecting  the  boundary  of  the 
spray  and  calculating  the  cone  angle  consists  of  the  following  steps: 

1 .  Capture  the  spray  using  a  high  speed  camera  and  save  the  avi  file  using  a  commercial  video  package. 

2.  Select  a  frame  from  the  recorded  avi  file  and  create  a  grayscale  intensity  image  as  well  as  zero  pad  the  signal. 

3.  Create  a  binary  image  from  the  grayscale  intensity  image  using  the  function  gray  thresh  using  a  commercial 
software  pacakge17.  This  function  utilizes  Otsu’s  method  which  chooses  the  threshold  to  minimize  the 
interclass  variance  of  the  black  and  white  pixels. 


4.  From  the  binary  image,  trace  the  spray  cone  boundary  where  non-zero  pixels  represent  the  boundary  of  the 
spray  cone  and  zero  pixels  constitute  the  white  background. 

5.  Detect  the  top  and  the  bottom  boundaries  of  the  spray  cone  at  a  desirable  downstream  location  from  the  exit 
nozzle. 

6.  Use  the  two  points  and  then  least-square  fit  a  line  parallel  to  the  boundaries  of  the  spray  cone,  and  find  an 
intersection  point  that  lies  on  the  axis  of  the  injector  to  calculate  the  spray  cone  angle.  In  order  to  determine 
the  axis  of  the  injector,  a  picture  was  taken  of  a  known  calibration  rod  that  was  placed  on  the  axis  of  the  exit 
nozzle. 

7.  Repeat  the  process  for  all  frames  for  the  total  duration  of  the  captured  video. 

8.  Plot  changes  in  total  cone  angle  as  a  function  of  time. 

9.  Repeat  the  process  for  different  downstream  locations  from  the  exit  nozzle  plane. 

Applying  the  steps  above,  oscillation  patterns  in  spray  cone  angles  as  a  function  of  time  can  be  determined. 

Similarly,  this  technique  is  used  to  measure  the  oscillation  of  the  air  core  diameter  as  a  function  of  time.  Note  that 

the  local  radial  distortions  are  irrelevant  since  the  experiment  focused  looking  at  a  planar  projection  of  the  cone. 

These  signals  are  then  analyzed  using  a  fast  Fourier  transform  of  the  spray  cone  angle  or  air  core  radius  histories10. 

Table  2:  High  speed  camera  settings  for  the  spray  and  air 
core  acquisition 


Parameter 


Nozzle  Exit 
Spray  Angle 


Vortex  Chamber 
Air  Core  Diameter 


Exposure  time  (microsec) 

Image  rate  (pps) 

ffl 

Duration  (sec) 

Image  width  (pix) 

Image  height  (pix) 

Image  resolution  (mm/pix) 
Number  of  images  per  test 


50 

63 

10000 

15037 

ffl  1.0 

ffl  1.0 

0.7207 

0.7188 

512 

512 

384 

256 

0.240 

0.126 

7207 

10810 

The  uncertainty  analysis  for  the  experiment  was  calculated  using  the  methodology  of  Coleman  and  Steele14.  For  the 
spray  cone  half  angle,  the  two  measured  variables  were  1)  the  downstream  distance  from  the  nozzle  exit  plane  and,  2) 
the  distance  measured  from  the  spray  boundary  to  injector  axis.  It  was  determined  that  the  uncertainty  for  spray  cone 
half  angle  between  40  and  50  degrees  was  less  than  2%.  Similarly,  uncertainty  for  typical  air  core  diameter 
measurement  ranging  from  0.25-0.5  inches  (6-12  mm)  was  less  than  6  %. 

III.  Results  and  Discussions 

The  swirl  injector  experiment  was  conducted  with  16  hole  dynamic  inlet  cap  at  three  ullage  pressures,  55,  70  and 
80  psia  (3.4,  4.8,  5.4  atm),  to  evaluate  the  effects  of  the  unsteady  flow  on  the  formation  of  the  spray  and  the  air  core 
diameter.  During  the  experiment  the  internal  flow  did  not  show  any  signs  of  two  phase  flow.  Pulsating  frequencies 
were  varied,  up  to  500  Hz,  by  adjusting  the  motor  speed.  For  all  ullage  pressures,  images  of  the  spray  cone  near  the 
exit  nozzle  and  the  air  core  inside  the  vortex  chamber,  as  well  as  the  pressure  readings  at  the  tangential  inlets  (pressure 
transducer  (k)  in  Fig.  2)  were  captured  and  recorded  simultaneously.  Figures  4  and  5  provide  typical  manifold  pressure 
measurements  for  the  long  and  short  injectors  respectively  at  a  tank  ullage  pressure  of  80  psi  (5.4  atm).  The  manifold 
sees  lower  pressure  then  the  ullage  pressure  because  of,  1)  the  pressure  drop  across  the  plumbing  system  and,  2)  the 
presence  of  the  dynamic  pressure  of  the  swirling  flow  created  by  the  dynamic  cap  in  the  manifold.  As  the  manifold  is 
of  limited  volume,  it  does  respond  to  the  transient  opening  and  closing  of  tangential  inlet  ports.  At  low  frequencies, 
the  manifold  pressure  oscillation  is  roughly  10-15%  of  the  mean  pressure,  while  at  higher  frequencies  the  oscillation 


amplitude  is  reduced  to  a  few  percent  of  the  mean.  There  do  not  appear  to  be  large  differences  in  the  magnitude  of 
these  oscillations  between  the  long  and  short  injectors. 

The  pressure  signals  were  analyzed  for  frequency  content  using  the  FFT  utility  in  a  commercial  software 
package17.  As  expected,  strong  peaks  were  found  at  the  low  driving  frequency  (61  Hz  in  long  injector,  50  Hz  in  short 
injector).  At  these  lower  frequencies,  the  second  harmonic  frequency  showed  the  next  highest  response  -  this  is 
evident  as  a  “beat  frequency”  in  both  of  the  signals  in  Fig.  4.  Some  activity  was  noted  at  the  subharmonic  tone,  but  at 
a  substantially  reduced  level. 

The  higher  frequency  cases  showed  an  interesting  dynamic  content.  For  both  the  160  Hz  case  with  the  long 
injector  and  140  Hz  case  with  the  short  injector,  the  subharmonic  (80  and  70  Hz,  respectively)  showed  FFT  response 
comparable  or  greater  than  that  of  the  primary  driving  frequency.  Negligible  response  was  noted  at  the  second 
harmonic  for  these  higher  frequency  cases,  the  short  waves  associated  with  this  harmonic  must  be  significantly 
damped  in  the  radial  passage  between  the  rotating  cap  and  the  manifold  outer  wall. 
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Figure  4:  Measured  manifold  pressures  (location  k  in  Fig.  2)  for  the  long  injector  at  two  different  frequencies, 

at  an  ullage  pressure  of  80  psi  (5.4  atm) 
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Figure  5:  Measured  manifold  pressures  (location  k  in  Fig.  2)  for  the  short  injector  at  two  different  frequencies, 

at  an  ullage  pressure  of  80  psi  (5.4  atm) 


The  average  mass  flow  produced  by  the  injector  was  determined  over  a  range  of  frequencies  using  a  standard 
“catch  and  weigh”  procedure.  For  all  three  ullage  pressure  settings  55,  70  and  80  psia  (3.7,  4.8,  5.4  atm),  the  motor 
was  set  to  four  different  speeds,  in  order  to  determine  the  mass  flow  rate  at  different  pulse  frequencies.  The 
experiment  was  run  and  the  spray  was  captured  for  60  seconds  to  calculate  the  mass  flow  rate  in  kilogram  per  second. 
The  resulting  measurement  has  a  maximum  uncertainty  of  1.7  percent.  Assuming  a  steady  flow  at  the  average 


manifold  pressure,  the  theoretical  mass  flow  rate  calculated  using  classical  swirl  injector  inviscid  theory6  was  0.476, 
0.558,  and  0.606  kg/sec  for  PuiiagQ  at  55,  70,  and  80  psia  (3.4,  4.8,  5.4  atm)  respectively  for  both  injectors. 

The  measurements  are  summarized  in  Figs.  6  and  7.  In  general  the  average  measured  mass  flow  rate  decreased 
slightly  as  frequency  increased,  presumably  due  to  the  dynamics  of  the  tangential  inlet  opening  and  closing  at  greater 
rates  as  inlet  cap  speed  increased.  Even  though  it  utilizes  the  same  inlet  and  nozzle  configuration,  the  long  injector 
generates  flow  rates  about  10%  higher  than  that  of  the  short  injector.  Qualitatively,  this  may  be  explained  by  the  fact 
that  longer  chamber  imposes  more  friction  on  the  rotating  flow,  hence  decreases  its  momentum,  which  makes  the  core 
diameter  to  decrease,  which  makes  the  cross  sectional  area  of  the  flow  in  the  nozzle  of  increase,  which  increases  the 
flow  rate19. The  discharge  coefficient  data  (Fig.  7)  show  values  above  unity  for  most  of  the  conditions  measured 
implying  that  measured  flows  are  higher  than  those  computed  on  a  1-D  inviscid  basis.  Since  the  theoretical  values 
should  represent  an  upper  bound,  it  is  speculated  that  there  may  have  been  some  water  hammer  effects  induced  in  the 
annular  passage  between  the  rotating  cap  and  the  outer  manifold  wall.  In  particular,  the  rotation  of  the  cap  imposes  a 
large  dynamic  response  at  the  inlet  to  the  tangential  channel  (a  main  factor  motivating  this  type  of  pulsator  design) 
which  would  presumably  be  substantially  larger  than  the  oscillation  detected  in  the  manifold.  The  postulated  water 
hammer  effect  is  stronger  on  the  long  injector  as  the  radial  passage  over  which  the  pulsations  occur  is  longer  in  extent 
as  well.  Manufacturing  issues  preclude  the  installation  of  a  pressure  tap  within  this  buried  location,  so  there  is  not 
sufficient  information  to  confirm  this  hypothesis,  however  the  manifold  pressure  measurements  (Figs  4,5)  do  show 
slightly  higher  pulsations  for  the  long  injector. 
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Figure  6:  Measured  average  mass  flow  rate  for  16  hole  dynamic  inlet  cap  for  the  short  and  the  long  swirl 

injector. 
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Figure  7:  Discharge  coefficient  behavior  for  the  short  and  the  long  swirl  injector. 


Spray  Angle  Measurements 

Over  a  broad  range  of  frequencies,  the  pulsations  within  the  conical  sheet  were  easily  viewed  with  the  naked  eye. 
Using  the  signal  processing  technique  outlined  in  the  prior  section,  quantitative  measurements  of  spray  cone  angle 
were  made  at  various  axial  stations.  Figure  8  shows  two  typical  images  captured  by  the  high  speed  camera  during  a 
forced  excitation  of  the  injector.  The  spray  cone  at  the  two  separate  times  shown  in  Fig.  8  has  cone  angle  varying  as 


much  as  7.9  degrees  when  measured  at  location  1.0  Dn  downstream  from  the  nozzle  exit  plane  (denoted  as  the  vertical 
line  in  the  images).  In  addition,  the  intensity  variations  of  the  light  passing  through  the  conical  films  indicate  that  the 
spray  film  thickness  also  changes  during  the  pulsed  state.  In  Fig.  8  the  captured  spray  seen  on  the  right  has  a  larger 
maximum  film  thickness  than  that  on  the  left.  The  high  amplitude  pulsations  at  this  modest  frequency  cause 
corrugations  in  the  sheet  as  evidenced  in  both  images. 


Figure  8:  Typical  images  of  the  conical  sheet  in  the  near  exit  region  of  the  orifice  using  the  short  injector  at  a 
pulse  frequency  100.1  Hz  and  Puiiage  ~  55  psia.  The  measured  spray  angle  is  94.6  deg  (a)  and  86.7  deg  (b). 


Figure  9  illustrates  the  spray  angle  variation  at  1.0  Dn  downstream  from  the  nozzle  exit  plane  on  the  left,  and  the 
corresponding  Fast  Fourier  Transform  (FFT)  signal  on  the  right.  The  short  injector  data  is  depicted  here  for  a  manifold 
pressure  of  55  psia.  For  frequencies  greater  than  30  Hz  ,the  spray  cone  angle  varies  in  a  saw-toothed  manner  as 
evidenced  in  the  waveforms  at  the  three  frequencies  shown.  If  the  total  angle  is  decomposed  into  two  halves  extending 
from  the  axis  of  the  nozzle  to  the  boundary  of  the  spray  cone,  the  top  and  the  bottom  half  of  spray  angle  generally 
show  to  be  in  phase. 

The  images  on  the  right  side  of  Fig.  9  are  the  FFTs  of  the  waveform  signals  on  the  left.  The  FFTs  show  a  strong 
peak  at  the  driving  frequency  whose  amplitude  is  dependent  on  the  driving  frequency  itself.  Large  amplitude 
pulsations  are  present  (well  within  the  accuracy  of  the  measurement)  over  the  frequency  range  shown.  There  is  little  if 
any  evidence  of  subharmonic  or  higher  harmonic  content,  although  the  latter  could  not  be  assessed  at  the  highest 
frequency  due  to  potential  aliasing  errors. 

Figure  10  shows  the  overall  cone  angle  response  curves  for  both  the  short  and  long  injectors  at  the  three  ullage 
pressures  used  in  the  study.  The  amplitude  of  the  cone  angle  fluctuation  is  not  the  same  for  all  driving  frequencies. 
Results  are  shown  using  two  distinct  measurement  stations  1  and  2.5  nozzle  diameters  downstream  of  the  orifice  exit 
plane.  Despite  the  existence  of  forces  such  as  surface  tension,  aerodynamic,  and  gravity,  the  major  portion  of  the  spray 
angle  pulsation  is  the  driving  force  from  the  pulsator.  In  general,  the  results  at  all  three  measurement  stations  show  the 
same  trend,  although  the  results  at  x=2.5  Dn  show  the  largest  amplitude  response  in  this  location  furthest  from  the  exit. 
The  maximum  response  occurs  near  300  Hz  for  the  short  injector  and  at  frequencies  near  250  Hz  for  the  longer 
injector.  The  difference  in  the  maximum  frequency  response  is  due  to  the  short  injector  having  a  higher  resonant 
frequency  than  the  long  injector.  Power  limitations  in  the  motor  driving  the  swirl  cap  did  not  allow  us  to  fully  capture 
the  peak  response  in  the  long  injector  at  the  highest  manifold  pressure.  In  general,  peak  response  frequencies  tend  to 
increase  with  ullage/manifold  pressure. 

There  is  also  some  evidence  of  a  peak  in  the  response  at  low  frequencies  below  100  Hz,  although  this  peak  is  less 
pronounced  than  the  main  resonance.  The  swirl  cap  is  rotating  at  fairly  low  speeds  at  the  lowest  frequencies  evaluated 
and  it  is  possible  that  the  pulsations  are  becoming  lower  in  amplitude  as  a  result.  Clearly,  the  cone  angle  response  (as 
measured  on  this  basis)  should  drop  to  zero  at  zero  frequency  as  we  would  expect  a  steady  conical  sheet  under  this 
condition.  For  this  reason,  the  measurement  is  not  completely  reflective  of  dynamic  response,  however,  it  is  believed 
to  be  a  good  indicator  at  the  higher  frequencies  where  waves  traveling  along  the  sheet  are  more  comparable  in  length 
to  the  sheet’s  thickness. 
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Figure  9:  Measured  spray  cone  angle  for  short  injector  at  1.0  Dn  for  Puiiage  ~  55  psia  (left).  Frequencies 

detected  from  applying  FFT  to  total  spray  cone  angle  data  (right).  Peak  amplitude  is  located  at  100.1  Hz  (top), 
300.3  Hz  (middle)  and  440.7  Hz  (bottom).  The  actual  pulsation  frequency  was  101.0  Hz  (top),  300.0  Hz  (middle) 

and  442  Hz  (bottom). 
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Figure  10:  Frequency  response  curves  for  summary  plot  of  spray  cone  magnitude  measured  at  1.0  and  2.5  Dn 
downstream  from  the  nozzle  exit  plane  at  Puiiage  55,  70  and  80  psia  for  the  short  injector  (left  column)  and  the 

long  injector  (right  column). 


In  order  to  understand  the  evolving  profile  of  the  spray  cone  generated  during  the  dynamic  state  a  series  of  images 
of  the  same  pulsed  flow  condition  were  captured  near  the  exit  nozzle.  Figure  11  shows  two  distinct  forms  of  spray 
being  generated.  The  image  on  the  left  is  shaped  like  a  trumpet  (convex),  with  the  spray  curing  outward  from  the  axis 
of  the  injector,  whereas  the  image  on  the  right  is  shaped  like  a  tulip  (concave)  with  the  spray  boundary  turning  inward 
to  the  axis  of  the  injector.  This  spray  angle  oscillation  is  due  to  changes  in  circumferential  velocity  inside  the  vortex 
chamber.  As  the  holes  in  the  dynamic  cap  rotate  about  the  holes  in  the  stationary  cap,  the  mutual  cross-sectional  area 
changes,  altering  the  velocity  of  the  flow  entering  tangential  inlets  (that  is  mass  flow  rate).  This  change  in  the  flow 
velocity  at  the  inlet  changes  the  circumferential  velocity  inside  the  vortex  chamber  and  propagates  downstream 
towards  the  nozzle  exit  and  results  in  the  oscillation  of  the  cone  angle.  Finally,  a  pulse  wave  is  observed,  where  the 
atomizing  jet  breaks  up  into  regularly  spaced  clumps  along  the  flow  direction20.  The  bunching  of  cones  and  waves 
seen  in  the  experiment  was  similar  to  the  findings  reported  in  previous  work15. 


Figure  11:  Gross  appearance  of  the  conical  sheet  boundary:  a)  Trumpet-shape  and  b)  tulip  conical  sheet  shapes 

at  P ullage  ~  70  psia  for  short  injector. 


Dynamic  Measurement  of  Air  Core  Diameter 


Another  set  of  experiments  was  performed  to  measure  the  air  core  diameter  fluctuation  inside  the  vortex  chamber 
during  the  dynamic  state.  Figure  12  shows  air  core  surfaces  at  two  different  times  while  the  pulsator  is  in  operation. 
There  is  a  clear  undulation  in  the  air  core  diameter,  visible  in  Fig.  12,  indicating  the  existence  of  a  wave  superimposed 
with  the  air  core  inside  the  vortex  chamber.  This  phenomenon  was  observed  for  frequencies  up  to  500  Hz,  which  was 
the  maximum  frequency  obtained  with  the  16  hole  dynamic  cap.  This  phenomenon  was  not  observed  in  a  non-pulsed 
system,  which  was  confirmed  by  the  cone  angle  fluctuation  diminishing  to  zero  in  Fig.  10 


Figure  12:  Air  core  inside  the  long  vortex  chamber  with  pulse  frequency  ~  484  Hz  and  Puuage  ~  55  psia.  Flow 

direction  is  left  to  right. 


The  pulsator  changed  the  vortex  chamber  air  core  diameter  as  much  as  4.5%  for  some  frequencies.  It  was  noted 
that  the  air  core  diameter  fluctuation  magnitude  maximized  at  certain  frequencies  indicating  a  resonance  similar  to  that 
observed  in  the  spray  angle.  Similar  image  processing  was  carried  out  on  the  air  core  diameter  measurements  and 
FFTs  of  these  signals  were  constructed  at  various  frequencies  to  assess  overall  response  characteristics.  These  results 
are  summarized  in  Fig.  13  for  both  short  and  long  injectors  at  ullage  pressures  of  55  and  70  psia  where  the  distance 
A350  is  at  an  axial  station  2.15  inches  (55  mm)  of  the  nozzle  exit  plane.  We  were  unable  to  make  measurements  at 
the  highest  ullage  pressure  due  to  difficulties  in  providing  adequate  power  from  the  motor  to  drive  the  swirl  cap  at  this 
highest  rotation  rate.  Figure  13  shows  resonant  conditions  at  frequencies  very  similar  to  those  detected  with  the  spray 
angle  measurements  (Fig.  9).  In  this  case,  there  is  no  indication  of  a  subharmonic  peak  and  the  amplitude  of  the 
pulsation  grows  toward  a  maximum  at  zero  frequency  as  one  would  expect  from  linear  theory.  There  is  some 
evidence  of  a  second  peak  in  the  55  psia  ullage  pressure  data  for  the  long  injector,  but  this  conclusion  is  tentative  as  it 
could  not  be  replicated  in  the  shorter  injector  and  because  there  were  only  a  few  data  points  establishing  this  trend.  In 
general,  the  trend  lines  in  Fig.  13  (and  Fig.  10)  are  more  for  clarity  in  capturing  the  data  than  in  suggesting  specific 
curvature/slope  in  a  given  region. 
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Figure  13:  Measured  dynamic  response  of  air  core  diameter  inside  the  vortex  chamber  at  Punage  55  and  70  psia 
for  the  short  injector  (left  column)  and  the  long  injector  (right  column). 


Comparison  of  Resonant  Frequencies  with  Theory 


Recently,  Ismailov  developed  an  analytic  expression  for  swirl  injector  resonance  based  on  the  wave  speeds  and 
volume  of  fluid  within  the  vortex  chamber.  The  methodology  follows  a  similar  path  as  that  used  to  find  resonant 
wavelengths/frequencies  for  waves  entering  a  deep  harbor  as  the  volume  within  the  harbor  interacts  with  waves 
entering  its  constricted  entry  much  like  waves  in  the  vortex  chamber  interact  with  the  nozzle  of  the  swirl  injector. 
Ismailov  and  Heister11  provide  a  derivation  of  this  simple  result,  referred  to  the  Abrupt  Contraction  Resonance  Model 
(ACRM).  The  conical  inlet  to  the  nozzle  is  necessarily  neglected  in  this  simple  result  and  the  vortex  chamber  is 
lengthened  by  lA  of  the  conical  entry  length  to  provide  an  approximate  dimension  for  wave  resonance.  The  predicted 
primary  resonance  frequency,  go,  from  this  model  is: 
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Where  Vt  is  the  injection  velocity  through  then  tangential  inlets,  Lv  is  the  effective  length  of  the  vortex  chamber 

including  half  of  the  nozzle  contraction  length,  and  rvc  is  internal  radius  to  the  surface  of  the  film  in  the  vortex 
chamber.  The  vortex  chamber  film  radius  and  tangential  inlet  velocity  are  coupled  via  the  pressure  drop  across  the 
film  as  presented  in  the  classical  swirl  injector  literature6,9.  This  result  indicates  that  shortening  the  vortex  chamber 


length,  Lvc  or  increasing  the  injector  flow  rate  (i.e.  Vt )  will  tend  to  increase  resonant  frequencies.  Using  this  result, 

resonant  frequencies  were  computed  for  each  of  the  three  ullage  pressures  for  both  long  and  short  injectors.  During 
this  process,  we  became  concerned  that  the  dynamic  swirl  cap  may  be  inducing  effects  not  considered  in  the  model  as 
evidenced  by  the  manifold  pressure  and  discharge  coefficient  data  in  Fig.  4-7.  There  appears  to  be  evidence  of  water 
hammer  effects  that  would  substantially  affect  the  inflow  pressure/velocity  as  a  result  of  the  periodic  opening  and 
closing  of  tangential  inlet  passages.  For  this  reason,  the  effective  feed  pressure  was  raised  in  an  effort  to  capture  this 
effect;  a  value  of  40  psia  seems  to  provide  an  excellent  result  as  indicated  in  Table  3.  In  fact,  the  correction  may  likely 
be  dependent  on  the  speed  of  the  rotating  cap,  but  there  is  insufficient  information  to  produce  a  logical  result  given  the 
manifold  pressure  is  the  only  measurement  and  this  is  far  displaced  from  the  tangential  inlet.  Clearly  this  is  an  area 
that  would  benefit  from  further  exploration. 

Table  3  provides  an  overall  comparison  of  the  experimental  measurements  (both  cone  angle  and  air  core  diameter) 
with  base  ACRM  and  corrected  ACRM  results.  The  experimental  measurements  are  in  good  agreement  (less  than 
5%)  between  the  two  approaches.  The  uncorrected  ACRM  values  are  significantly  lower  than  those  in  the 
experiments.  Applying  the  40  psia  correction  brings  the  ACRM  values  very  close  to  the  measured  results.  While  we 
believe  that  there  is  some  basis  for  the  correction,  clearly  more  measurements  are  required  to  fully  assess  the  dynamic 
character  of  the  swirl  cap  and  its  influence  on  the  massflow  pulsation  produced  at  various  rotation  rates. 


Table  3:  Comparison  of  primary  peak  response  frequencies  (in  Hz)  of  measured  nozzle 
exit  spray  angle  and  values  computed  from  the  ACRM 


Short  Injector 

Long  Injector 

P  ullage 

55  psia 

70  psia 

80  psia 

55  psia 

70  psia 

80  psia 

Measured  Peak  Spray 
Angle 

Measured  Peak  Vortex 

305 

324 

345 

241 

261 

n/a 

288 

307 

— 

242 

274 

Core  Diameter 

ACRM,  uncorrected 

205.7 

240.9 

261.8 

170.7 

200.0 

217.3 

ACRM,  corrected 

290.3 

316.3 

332.4 

241.0 

262.6 

276.0 

Generation  of  Travelling  Waves 

A  series  of  images  of  a  propagating  wave  inside  the  long  injector  vortex  chamber  for  excitation  frequency  set  at 
491  Hz  is  shown  in  Fig.  14.  Note,  the  two  arrows  on  all  images  represent  the  location  of  the  maximum  air  core 
diameter  in  each  frame. 


Time  =  0.201 97  sec 


Time  =  0.2023  sec 


Time  =  0.20263  sec 


Figure  14:  Series  of  images  of  a  propagating  wave  inside  the  long  injector  vortex  chamber  for  Puiiage  55  psia 
with  pulse  frequency  set  at  491  Hz.  Flow  direction  is  left  to  right. 


The  following  phenomena  are  observed.  In  frame  (a),  a  surface  wave  is  located  at  the  head  end  of  the  vortex 
chamber.  A  surface  wave  which  is  generated  by  a  flow  disturbance  at  the  head  end  propagates  downstream  along  the 
injector  axis,  as  captured  in  frames  (b)  through  (e).  The  wave  is  travelling  downstream  and  into  the  converging  part  of 
the  vortex  chamber  while  rotating  about  the  axis  of  the  air  core  with  some  large  amplitude  as  captured  in  frames  (g) 
and  (h).  When  the  wave  reaches  the  nozzle,  part  of  it  goes  into  the  nozzle,  part  of  it  reflects  back.  The  wave  that  exists 
in  the  vortex  chamber  is  a  result  of  superposition  of  forward  and  backward  traveling  waves.  Therefore,  the  observed 
wave  length  in  the  vortex  chamber  is  not  exactly  equal  to  the  wavelength  that  corresponds  to  491  Hz.  A  newly 
generated  wave  then  starts  again  from  the  vortex  chamber  head  end  and  repeats  the  process.  The  wave  images  showed 
that  the  disturbance  wave  characteristics  depend  on  the  disturbance  frequency,  as  this  effect  is  occurring  at  the 
excitation  frequency  of  the  pulsator.  The  period  of  this  travelling  wave  approximately  matches  the  pulsation 
frequency.  That  is,  1/(0.20463-0.2023)  =  430  Hz.,  and  pulse  frequency  is  at  491  Hz. 

IV.  Conclusions 

An  experimental  study  of  the  dynamics  of  a  classical  (simplex)  swirl  injector  has  revealed  the  behavior  of  the 
conical  sheet  formed  near  the  injector  exit  and  internal  film  within  the  vortex  chamber.  A  super-sized,  transparent 
injector  configuration  was  used  with  a  unique  pulsator  design  that  periodically  blocked  and  opened  tangential  inlet 
channels  and  in  doing  so  provided  high  amplitude  perturbations.  A  16  hole  rotating  cap  was  used  to  create  the 
periodic  blockage  and  measurements  were  successfully  taken  at  frequencies  as  high  as  500  Hz.  Manifold  pressure 
measurements  do  indicate  unsteadiness  due  to  the  pulsations  induced  by  the  operation  of  the  rotating  cap  and  the  mean 
discharge  characteristics  of  the  device  reveal  that  there  may  be  some  water  hammer  induced  flow  as  the  discharge 
coefficients  exceed  unity  under  most  operating  conditions. 

Observation  of  the  dynamics  of  the  conical  sheet  in  the  region  near  the  nozzle  exit  provided  an  excellent  means  of 
assessing  the  dynamic  response  of  the  injector.  Data  were  repeatable  at  various  axial  stations  and  provided  good 
signal/noise  characteristics  over  a  broad  range  of  frequencies.  Resonant  behavior  was  observed  in  these  data  as 
theorized  in  the  recent  work  of  Ismailov9.  Comparisons  of  the  measured  resonant  frequencies  with  the  theoretical 
results  of  Ismailov  and  Heister11  yielded  poor  results,  but  if  the  theory  used  an  adjusted  manifold  pressure  to  reflect 
potential  water  hammer  effects  an  excellent  correlation  was  obtained.  Unfortunately,  the  experimental  apparatus  was 
not  amenable  to  further  exploration  of  water  hammer  effects  and  this  would  be  a  topic  for  further  study. 

Finally,  traveling  waves  of  the  type  theorized  by  Bazarov6  were  imaged  within  the  vortex  chamber  at  non-resonant 
conditions.  The  speed  of  these  waves  correlated  fairly  well  with  the  pulsation  frequency. 
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